A Mathieu function boundary spectral method for scattering by multiple variable poro-elastic plates, with applications to metamaterials and acoustics

1. Introduction

Motivated by many applications, there is substantial interest in solving Helmholtz scattering problems on unbounded domains with complicated boundary conditions. In this article, we consider the situation of Helmholtz scattering off (multiple) finite plates in two dimensions. When embedded in three dimensions, this corresponds to plates of infinite span but finite chord. When the geometry and boundary conditions are sufficiently simple, a successful approach for this problem is the Wiener–Hopf method [1–3]. For example, the Wiener–Hopf method allows one to capture the interaction of a semi-infinite edge with a quadrupole source and compute the far field. However, typically in such situations, one would want to model the interaction between the leading and trailing edges of a finite plate, which is important as both backscattering of the trailing-edge field by the leading edge [4] and structural resonances can be significant. There are some extensions of the Wiener–Hopf method which can deal with finite plates, but such extensions are non-generic and difficult due to the need to solve a matrix, rather than a scalar, Wiener–Hopf equation. Another common case encountered in applications, which cannot be tackled by the Wiener–Hopf method, is when physical parameters vary along the boundary of the domain. Such variations are expected to be crucial in biological applications [5] and to avoid discontinuous boundary conditions where additional scattering occurs [6]. Variation in physical parameters is also important in the study of metamaterials, such as acoustic black holes (see §5), which rely on a smooth variation of stiffness that, in the right circumstances, leads to almost 100% absorption of the incident wave energy [7,8]. Interactions of acoustic or hydrodynamic fluctuations with thin elastic structures arise in numerous other situations such as aerodynamic noise reduction [6,9–11] and the modelling of ice sheets and marine platforms in oceanography [12–15]. In all such cases, accurate and fast numerical methods are key to predicting the effect of external forces and variable parameters such as elasticity on an elastic plate, or the effect of elasticity on the radiated field, and thus crucial for providing insight into a wide range of fluid dynamical problems.

By starting with separation of variables in elliptic coordinates, we develop a boundary spectral method for scattering by multiple variable poro-elastic plates. This allows both accurate and rapid computation of the scattered field, as well as great flexibility in the boundary conditions specified on the plates. Separation of variables leads to angular Mathieu equations and radial Mathieu equations, and the solutions to these equations are the well-known Mathieu functions [16,17]. Historically, the problem of plane wave scattering of a rigid screen was first rigorously studied by Schwarzschild [18] based on the Sommerfeld half-plane problem and shortly after by Sieger [19] by employing Mathieu functions. Some numerical work based on this solution was presented in [20,21], and more recently in [22]. Extensions with different boundary conditions on elliptic shells were considered in [22,23]. Mathieu functions were also shown to be an effective tool for low-frequency scattering of a rigid (non-porous) plate in [24], where comparisons were made with semi-analytical boundary integral methods.

This article demonstrates that Mathieu functions offer a direct and rapid approach to tackle many interesting boundary value problems. To the authors’ best knowledge, the problem of acoustic scattering from multiple elastic plates with varying elasticity (or even a single plate with varying elasticity) using Mathieu functions has not been treated before. Our solution representation directly provides a sine series approximation of the far-field directivity and, unlike standard boundary methods, is easy to evaluate near the scatterers. This means that the near field can be computed efficiently and in a stable manner. These advantages mean that it is particularly good for a simple model of turbulence using Lighthill’s analogy. For example, the numerical method allows rapid and easy calculation of structural or acoustic resonances, which are generally challenging to compute when the physical parameters vary along the plate [25–27] or when sophisticated plate theories are involved [28]. To demonstrate the flexibility of the local Mathieu function expansions for arbitrarily positioned plates in two dimensions, figure 1 shows the total field for a quadrupole source scattering off four elastic plates. We also note that boundary conditions additional to those considered in §2 can easily be incorporated. It is important to point out, however, that the approach of this article cannot deal with curved boundaries which do not have a local coordinate system in which to perform separation of variables. Code for the numerical method is provided at https://github.com/MColbrook/MathieuFunctionCollocation.

• Download figure
• Open in new tab
• Download PowerPoint

Problems similar to a poro-elastic finite plate include the case of semi-infinite plates that are uniformly porous [1], or uniformly poro-elastic [11], which can be treated using scalar Wiener–Hopf techniques. These examples can be extended to more complicated porous boundary conditions [6,29], but in such cases, the analysis leads to a matrix Wiener–Hopf equation which is more difficult to solve. Elastic properties have also complexified previous numerical simulations. For example, recent work [30] (extended to three dimensions in [31]) for the scattering of a near-field source by a finite perforated elastic flat plate requires two problems to be solved; one for the structural modes of the plate which is done via a spectral method; the second for the scattering of the acoustic source which is achieved via a boundary element method (BEM). We compare our results (in the restricted case of constant porosity and constant elasticity dealt with in [30]) to those of [30] in §3c, demonstrating that separation of variables yields a faster, more robust and more accurate method for the case of a single plate. See also [32,33] for an expansion scheme of the plate deformation connected to Chebyshev polynomials that tackles the problem of a single elastic plate in a rigid baffle (our numerical scheme can handle this problem with an appropriate modification of the boundary conditions when we separate variables in §3). Another approach for these types of problems is the unified transform [34] (see also [35–38] for recent developments), a Fourier space boundary spectral collocation method which in certain cases generalizes the Wiener–Hopf method [9,39].¹ However, using the unified transform in unbounded domains requires the setting up of several global relations by hand, which becomes complicated in complex geometries. More broadly, there has been recent interest in spectral methods to solve scattering problems that can be recast as a Riemann–Hilbert problem [43,44], though, as far as the authors are aware, such methods have not yet been applied to elastic or porous scatterers.

The structure of this article is as follows. In §2, we describe the mathematical model for a single plate. The numerical method is presented in §3, where we also compare with the boundary element method of [30]. Examples of diffraction by elastic plates of varying stiffness are presented in §4, including the peculiar effects of an acoustic black hole in §5; we are not aware of any previous studies of this effect in such plates [8]. In §6, we describe how to extend the method to multiple plates. Concluding remarks are given in §7.

2. Mathematical model for single plate

Suppose that an incident sound wave travels towards a plate situated at −d ≤ x ≤ d (where d > 0) and y = 0. The incident field will be denoted ϕ_I and the scattered field by ϕ. The incident pressure field is given by $p_{\text{I}}={\rho }_f{c}_0^2{\varphi }_{\text{I}}$, where ρ_f is the mean fluid density and c₀ the speed of sound, so that throughout we deal with dimensionless fields ϕ_I and ϕ. We assume that ϕ has the usual time dependence e${\hspace{0.17em}}^{-\text{i}\omega t}$ (omitted throughout) and therefore satisfies the Helmholtz equation

We consider poro-elastic boundary conditions. Other types of boundary conditions can also be tackled by the methods of this article (see, for example, the list of boundary conditions and physical interpretations in [45]), including non-local boundary conditions, but we stick to the following case for brevity. For completeness, we have also provided an electronic supplementary material detailing the implementation for rigid porous plates.

We consider a poro-elastic plate with evenly-spaced circular apertures of radius R, Rayleigh conductivity of K_R = 2R, and fractional open area α_H = NπR² (where N is the number of apertures per unit area) [46]. The plate deformation is given by $\eta (x){\text{e}}^{-\text{i}\omega t}$ (the time dependence is again assumed and omitted) and η(x) satisfies the thin-plate equation

3. Single plate solution

(a) Expansion of solution in Mathieu functions

The solution ϕ is an odd function in the variable y and hence we can consider solving the PDE system in the upper-half plane {(x, y):y > 0}. First, we introduce elliptic coordinates via x = dcosh (ν)cos (τ), y = dsinh (ν)sin (τ), where, with an abuse of notation, we write functions of (x, y) also as functions of (ν, τ). Elliptic coordinates for d = 1 are displayed in figure 2. The appropriate domain then becomes ν ≥ 0 and τ ∈ [0, π]. To simplify the formulae, we let $Q=d^2{k}_0^2/4$. Separation of variables leads to the expansion

$$\varphi (\nu ,\tau )=\sum _{m=1}^{\mathrm{\infty }}{a}_m{\text{se}}_m(\tau ){\text{Hse}}_m(\nu ),$$ 3.1

where se_m(τ) = se_m(Q;τ) denote sine-elliptic functions and Hse_m(Q;ν) = Hse_m(ν) denote Mathieu–Hankel functions. A full derivation is provided in the electronic supplementary material.

• Download figure
• Open in new tab
• Download PowerPoint

The functions se_m are expanded in a sine series as

$${\text{se}}_m(\tau )=\sum _{l=1}^{\mathrm{\infty }}{B}_l^{(m)}\sin (l\tau ).$$ 3.2

This Fourier series converges absolutely and uniformly on all compact sets of the complex plane [17] and we find the coefficients $B_l^{(m)}$ via a simple Galerkin method. The convergence to the eigenvalues and eigenfunctions depends on the parameter Q, in general being slower for larger Q. However, the convergence is exponential, yielding machine precision for small truncation parameter n, even for very large Q [42].

The functions Hse_m(ν) can be expanded using Bessel functions [16,17]:

$$\begin{array}{rl}{\text{Hse}}_m(\nu )& =\sum _{l=1}^{\mathrm{\infty }}\frac{{(-1)}^{l+m}{B}_l^{(m)}}{C_m}\\ & \times [J_{l-1}({\text{e}}^{-\nu }\sqrt{Q})H_{l+p_m}^{(1)}({\text{e}}^{\nu }\sqrt{Q})-J_{l+p_m}({\text{e}}^{-\nu }\sqrt{Q})H_{l-1}^{(1)}({\text{e}}^{\nu }\sqrt{Q})],\end{array}$$ 3.3

where p_m = 1 if m is even and p_m = 0 if m is odd. Here J_n denotes the Bessel function of the first kind of order n and we remind the reader that $H_n^{(1)}$ denotes the Hankel function of the first kind of order n. The series in (3.3) converges absolutely and uniformly on all compact sets of the complex plane [17]. We choose the normalization constants C_m such that Hse_m′(0) = 1. The terms in the series (3.3) can easily be evaluated for small l. However, for large l, the terms in the series suffer from underflow and overflow associated with cancellations between the Bessel and Hankel functions. For large l and fixed $x\in {\mathbb{R}}_{>0}$, we use the asymptotics

$$J_l(x)=\sum _{j=0}^q\frac{{(-1)}^j}{j!(j+l)!}{\left(\frac{x}{2}\right)}^{2j+l}+\mathcal{O}\left(\frac{1}{(q+l+1)!}\right)$$

and

$$H_l^{(1)}(x)=\frac{-\text{i}}{\pi }{\left(\frac{2}{x}\right)}^l\sum _{j=0}^q\frac{(l-j-1)!}{j!}{\left(\frac{x}{2}\right)}^{2j}+\mathcal{O}((l-(q+2))!),$$

valid as l → ∞. For fixed $a,b\in \mathbb{Z}$, this gives the asymptotic form

$$\begin{array}{rl}& J_{l+a}({\text{e}}^{-\nu }\sqrt{Q})H_{l+b}^{(1)}({\text{e}}^{\nu }\sqrt{Q})\\ & =\frac{-\text{i}}{\pi }{\left(\frac{\sqrt{Q}}{2}\right)}^{a-b}{\text{e}}^{-\nu (2l+a+b)}\left[\sum _{j=0}^q\frac{{(-1)}^j(l+a)!}{j!(j+l+a)!}{\left(\frac{{\text{e}}^{-\nu }\sqrt{Q}}{2}\right)}^{2j}\right]\\ & \times \left[\sum _{j=0}^q\frac{(l+b-j-1)!}{j!(l+a)!}{\left(\frac{{\text{e}}^{\nu }\sqrt{Q}}{2}\right)}^{2j}\right]+\mathcal{O}(l^{-(q+2)})\end{array}$$

We found this to be an excellent approximation for large l. It can also be accurately evaluated for moderate q since the terms (l + a)!/(j + l + a)! and (l + b − j − 1)!/(l + a)! can be evaluated as products of j and |j + 1 + a − b| terms, respectively. In what follows, we typically used this asymptotic form when l > 100 and took up to q = 5 terms. When plotting errors of our method, we were careful to compare against converged computations for which the series (3.3) was evaluated directly using extended precision (such checks were the only place where we made use of extended precision). Figure 3 shows the first 10 eigenfunctions and Mathieu–Hankel functions for k₀ = 20 computed to machine precision.

• Download figure
• Open in new tab
• Download PowerPoint

We use the boundary conditions to solve for the unknown coefficients a_m, after which the solution can be evaluated anywhere in the (x, y) plane. Of particular interest is the far-field directivity, D(θ), which is defined via

$$\varphi (r,\theta )\sim D(\theta )\frac{{\text{e}}^{\text{i}wr}}{\sqrt{r}},\text{as }r\to \mathrm{\infty },$$ 3.4

where (r, θ) are the usual polar coordinates. Given the Bessel function expansion of Hse_m(ν) in (3.3), we can directly compute D(θ) from (3.1) using asymptotics of Bessel functions (for large arguments, not large order as was used previously above). In the appropriate limit, τ becomes the polar angle θ, whereas ν becomes cosh ⁻¹(r/d) (in the far field the confocal ellipses can be approximated by concentric circles, see figure 2). We therefore have

$${\text{Hse}}_m(\nu )\sim \frac{{(-1)}^{m+1}}{C_m\sqrt{\pi r\sqrt{Q/d^2}}}\exp \left([2r\sqrt{Q/d^2}-\frac{(p_m+1)\pi }{2}-\frac{\pi }{4}]\text{i}\right)B_1^{(m)},\text{as }r\to \mathrm{\infty },$$ 3.5

and hence

$$D(\theta )=\sqrt{\frac{2}{\pi k_0}}\sum _{m=1}^{\mathrm{\infty }}\frac{a_m{B}_1^{(m)}}{C_m}\exp \left(\frac{(2p_m-3)\pi }{4}\text{i}\right){\text{se}}_m(\theta ).$$ 3.6

An advantage of our approach is that, for the case of a single plate, we implicitly compute a sine series for the far-field directivity D(θ) through the sine-elliptic functions se_m(θ) given by (3.2).

(c) Comparison with elastic boundary element method

In this section, we analyse the numerical performance of the proposed method with constant physical parameters. Further examples where parameters vary will be given in later examples. A comparison between our method and the unified transform for a rigid porous plate can be found in [42].

We compare the proposed collocation method with the BEM of [30], which deals with constant porosity and elasticity. The method of [30] first computes the spectral modes of the fourth-order derivative operator (acting on the left-hand side of (2.1)), before recasting the boundary conditions in terms of these vibration modes of the plate, and then solving the resulting boundary element scheme. In this section, we shall be consistent with the set-up of [30] and consider a plate that lies along {(x, 0):x ∈ [0, 1]}, is clamped at x = 0, and free to move at x = 1. To compare with the parameters of [30], for a plate of mass m per unit area and effective plate stiffness $\overline{B}$, we define² the coincidence frequency

$${\omega }_c={\left(\frac{(1-{\alpha }_H)m{c}_0^4}{\overline{B}}\right)}^{1/2},$$

the vacuum bending wave Mach number

$$\mathit{\Omega }={\left(\frac{\omega }{{\omega }_c}\right)}^{1/2}=\frac{k_0}{k_B},$$

and the intrinsic fluid-loading parameter

$$ϵ=\frac{{\rho }_f{k}_0}{(1-{\alpha }_H)m{k}_B^2}.$$

Note that since we are considering constant parameters in this subsection, ω_c, Ω and ϵ are constant. After a suitable rescaling with the plate length (and, with an abuse of notation, keeping the same notation for physical parameters), the non-dimensionalized boundary conditions become

$$(1-{\alpha }_H)\frac{{\mathrm{\partial }}^4\eta }{\mathrm{\partial }{x}^4}-\frac{k_0^4}{{\mathit{\Omega }}^4}\eta =-(1+\frac{4{\alpha }_H}{\pi })\frac{ϵ}{{\mathit{\Omega }}^6}{k}_0^3[\varphi ]$$ 3.9

and

$$\frac{\mathrm{\partial }\varphi }{\mathrm{\partial }y}{|}_{y=0}+\frac{\mathrm{\partial }{\varphi }_{\text{I}}}{\mathrm{\partial }y}{|}_{y=0}=(1-{\alpha }_H)k_0^2\eta +\frac{2{\alpha }_H}{\pi R}[\varphi ].$$ 3.10

A broad parametric study of how our collocation approach compares to BEM would be an exhaustive task. Instead, we provide some comparisons pertinent to the general performance of both methods. We therefore set R = 10⁻³, α_H = 2 × 10⁻³ and ϵ = 0.0021 throughout this section (representative of an aluminium plate in air [46]). We compare both methods for computing the far-field directivity (computed by measuring the scattered field at radius 100 for BEM), using a discrete relative L² error defined by

$$\sqrt{\frac{\sum _i|\tilde{D}({\theta }_i)-D({\theta }_i){|}^2}{\sum _i|D({\theta }_i){|}^2}}.$$

Here, $\tilde{D}$ is the computed directivity and D the true directivity which is estimated via a converged computation with larger M and N for our method, and a larger number of boundary elements and modes for BEM. The θ_i are taken to be 201 equally spaced points covering the interval [0, π]. Consistent with [30], we consider the case of placing a quadrupole at (x, y) = (1, 0.01) and compute the resulting far field of ϕ. We chose to compare the accuracy of computing the far-field directivity as opposed to the jump in pressure across the plate since the numerical approach of [30] adopts a small but positive plate thickness (however, we also obtain similar qualitative results for other physical quantities of interest). Therefore, we do not expect exact agreement between the BEM and our collocation approach (which deals with plates of zero thickness).

Figure 4 shows |D(θ)| for various Ω and k₀. These show excellent agreement between both methods (we used M = N for the Mathieu function collocation method). There is a slight deviation for k₀ = 20 and Ω = 0.05 due to the non-zero plate thickness in BEM (this is expected to make more of a difference for larger k₀ and smaller Ω). Figure 5 shows the convergence of BEM (default 100 modes) as a function of the number of degrees of freedom of the linear system. We see quite slow algebraic convergence (typical of standard BEM). For small k₀, the errors are smaller for larger Ω as the plate becomes more rigid. This was less pronounced for larger k₀. However, in this case, for smaller Ω we needed a larger number of modes for the error not to plateau. This is expected since, as a rough heuristic, the number of modes needed scales as the bending wavenumber k_B = k₀/Ω. Figure 6 shows the convergence of our Mathieu function collocation method, where we have also plotted the bending wavenumbers. For each set of parameters, there is a region of algebraic convergence (roughly cubic) once the number of degrees of freedom is of the order k_B. There is also an initial region of rapid convergence (typical of spectral methods) most pronounced for larger Ω. The Mathieu function approach achieves errors several orders of magnitude smaller than BEM and for much fewer degrees of freedom.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

Finally, figure 7 shows the average times of the methods implemented on a 5-year-old laptop, including to evaluate the far field. The Mathieu function approach is much faster (see the different scales on the vertical and horizontal axes), even when the size of the linear systems are the same. A possible reason for this is the implementation of the BEM code, however, as demonstrated in figures 5 and 6, much smaller system sizes are needed for a given accuracy when using the collocation method. For BEM, we have shown separately the times taken to compute the vibrational modes and also to set up and solve the linear system. When using BEM, the vibrational modes do not need to be recomputed for different parameters (assuming enough modes are included to capture the oscillations). However, the precomputation of the coefficients in the expansion (3.2) via a symmetric tridiagonal eigenvalue, which needs to be performed for each value of k₀ in the Mathieu function approach, takes negligible time compared to solving the linear system for large M.

• Download figure
• Open in new tab
• Download PowerPoint

Though we do not repeat the results here, a relative accuracy of approximately three digits and runtime of a few seconds (including evaluation) was reported in [9] for similar parameters using the unified transform. Therefore, our approach in this article is also faster and more accurate than the unified transform implemented in [9,49].

4. Diffraction by an elastic plate of varying thickness

For the rest of this article, we consider the choices

Here, we investigate how different variations in the plate thickness h(x) influence the scattered field. We define a functional P, proportional to the total above-plate scattered sound power

(a) Linear variation

Consider first a linear variation in the plate thickness for a plate clamped at both endpoints with

$$h(x)=0.004(1-cx),$$ 4.4

for different c. The sound power P is shown in figure 8a for an incident plane wave of angle π/3. Looking at the constant stiffness (c = 0), the most apparent feature is the presence of resonance peaks. These resonant contributions have decreasing power and increasing frequency width as the frequency increases [51]. The resonance peaks of a fluid-loaded plate are known to be slightly lower than the in vacuo plate resonance wavenumbers [30]. We have found that varying c does not significantly change the sound power P variation for different acoustic wavenumbers. The only observed effect is the shift in the resonance peaks, which is small for small c and larger for larger c (as expected, smaller modifications of h lead to smaller changes in P). Hence it is found that, similarly to the introduction of fluid-loading, linear variation in h (and monotonic changes in B) changes the position of the resonance peaks in figure 8.

• Download figure
• Open in new tab
• Download PowerPoint

5. Acoustic black hole

We next consider the case of an acoustic black hole. These are new physical objects, introduced and investigated over the last 15 years or so [7,8,52–56], that under certain circumstances can absorb almost 100% of the incident wave energy. Acoustic black holes have been investigated mainly for flexural waves in thin plates, where the local thickness varies according to a power law, with the power-law exponent being greater than or equal to 2. Here, we explore their properties in acoustic scattering. Whereas previous work considers incident waves that originate inside the plate/wedge, we consider the interaction of such a plate with an incident field.

(a) Incident plane wave

In this example, we take α_H ≡ 0 (i.e. zero porosity), and take the plate to be clamped at both endpoints. The thickness is chosen to vary according to

$$h(x)=0.001x^2+h_0,$$ 5.1

for a small positive cut-off h₀. If h₀ = 0, the thin-plate equation (2.1) becomes singular at x = 0. For this reason, and to also avoid physically impractical cases, we consider examples of small but positive non-zero h₀. We consider an incoming incident plane wave of angle 3π/4 and k₀ = 20. Figure 9 shows the plate deformations for h₀ = 10⁻⁶ and h₀ = 10⁻³. For small h₀, the plate vibrations become very large as the thickness decreases at x = 0. The oscillations become clustered near the thin portion of the plate (see the magnified section), and this effect is increased by making h₀ smaller. This effect is removed when h₀ = 10⁻³. Figure 10 shows the corresponding near fields. We see that near x = 0, the incident field is able to pass through the plate, causing little reflection for h₀ = 10⁻⁶. Again this effect is removed for the larger h₀ = 10⁻³. In figure 11, the far field is presented. There is a slight reduction in the scattered noise for smaller h₀, with a less focused scattering direction. Finally, figure 12a plots the convergence of the physical variables of interest and demonstrates that we can easily gain several digits of relative accuracy, even for small h₀.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

(b) Quadrupole sound source

The noise generated by the turbulence at the trailing edge of an aerofoil can make a significant contribution to the overall production of aeroacoustic noise, especially at high frequencies [57,58]. By Lighthill’s analogy, turbulent eddies are represented by a distribution of quadrupole sources in the same volume [59]. This motivated the study of a simplified model of the scattering by a plate with forcing given by a quadrupole at (x, y) = ( − 1, 0.001). The resulting near and far fields can be used to study aerofoil edge adaptions [30].

In this example, we take α_H ≡ 0 (i.e. zero porosity), and take the plate to be clamped at x = 1 but free at x = −1. The thickness is chosen to vary according to

$$h(x)=0.001{(x+1)}^2+h_0,$$ 5.2

for a small positive cut-off h₀. We consider the case of k₀ = 25 for h₀ = 10⁻⁶ and h₀ = 10⁻³. Convergence of the method is shown in figure 12b. We also checked our results near the vibrating tip by resorting to reciprocity. The value of the fields at (x, y) = ( − 1, 0.001) with a quadruple at (1, 1) is the same as the value at (x, y) = (1, 1) with the quadruple at ( − 1, 0.001).

Figures 13–15 show the plate deformations, near field and far field, respectively. The plate deformations behave qualitatively as before, with oscillations clustering near the thin part of the plate for smaller h₀. The imaginary part of η for h₀ = 10⁻⁶ is not zero, but it is small in comparison with the real part. This is because the real part of the incident quadrupole dominates near the source. We see a very interesting effect for the near field. The magnitude of the field is much smaller for h₀ = 10⁻⁶, and in fact appears to be dominated locally around the right tip (1, 0) which is unusual. We can also see that there are evanescent pressure waves on the surface of the plate in figure 14 (magnified). The more flexible end of the plate absorbs (rather than scatters) the pressure fluctuations and propagates them down the plate to the less flexible endpoint. On reaching the x = 1 tip, the pressure fluctuations scatter resulting in a directivity pattern as if the main source was located near the x = 1 endpoint. By contrast, for h₀ = 10⁻³, expected cardioid directivity around the point ( − 1, 0) is observed typical for such problems, and no evanescent pressure waves are visible. The corresponding pattern is observed in the far-field directivity, where we see that for h₀ = 10⁻⁶, the scattered field is reflected back in the direction of the source and is much smaller than that of h₀ = 10⁻³. The scattered field in the direction of the incident field is an interesting example of an acoustic black hole effect in a plate of varying elasticity. The authors are not aware of this effect being studied in such plates. The usual setting for this is an elastic wedge, where the cross-sectional thickness varies according to a power law [8]. One interesting observation is that the black hole effect relies on the power function going to nearly zero [8]. It is mitigated when the cross-sectional thickness decreases to 10⁻³ in figure 15b, where the majority of scattering obeys the usual reflection.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

The above results indicate that there is potential to exploit the acoustic black hole effect in edge adaptations. For example, acoustic black holes can be used to direct the vibrations away from the aerofoil edges towards the middle of the aerofoil where vibration-absorbing mechanisms can be placed. The need for small h₀ gives practical limitations for its use, but the truncated profiles with correctly placed damping layers can still be practical. There are some preliminary experimental studies for sound absorption in air [8]. However, currently, there is very little known about the use of acoustic black holes in aeroacoustics, and further theoretical and experimental investigations and validation are needed.

6. Extension to multiple plates

We can also use the numerical method in §3 to compute the scattered field from multiple plates. Suppose that we have plates P_[i] for i = 1, …, S, whose lengths are 2d_[i]. We also suppose that the open set ${\mathbb{R}}^2\mathrm{\setminus }({\cup }_{i=1}^S{P}_{[i]})$ is connected (in particular, we exclude the possibility that plates enclose a region, though this can be dealt with via suitable modifications). We use sub/superscripts [i] to denote quantities associated with the plate P_[i]. Each plate P_[i] induces a corresponding scattered field given by

Numerically, we solve this problem in the same way, where we take M_[i] Mathieu functions for the expansion along the ith plate and we supplement the expansion of ϕ_[i] with an expansion of η_[i] in terms of N_[i] Chebyshev polynomials of the first kind along the plate P_[i]. The relation (2.1) becomes

As a simple example, we consider the case of two plates where P_[1] is elastic and clamped with endpoints ( ± 1, L) and P_[2] is rigid with endpoints ( ± 1, − L). For P_[1], we set

Figure 16 shows the far fields for L = 0.01 and L = 0.5. In the acoustic compact case $L\ll k_0^{-1}$, the scattered field behaves as if it is incident on a single plate. This gives a symmetric scattered field, which does not vary monotonically with B (as expected due to effects such as resonances). For larger spacings, each edge (four in total) scatters an acoustic field which interacts in the far field to create an oscillatory directivity pattern. If the elastic plate is suitably flexible to be excited by the incident wave and absorb energy, its scattering will be distinctly different to a rigid plate, and hence alter the overall far field directivity. The primary effects are noticeable in the Fresnel lobes. Figure 17 shows the near field and far field for $L=0.1=\mathcal{O}(k_0^{-1})$. In this case, the plates support a specific ‘duct’ mode between them. The scattering of these modes by the edges contributes to the far-field noise. Altering the elasticity of the upper plate alters the fundamental structure of the duct and what modes can exist there. This too impacts the scattering in addition to direct scattering by each of the four edges.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

7. Conclusion

This article developed a boundary spectral method, based on collocation of local Mathieu function expansions, for Helmholtz scattering off multiple variable poro-elastic plates. Such boundary conditions are challenging for current methods, and we compared our approach to an elastic boundary element method in §3c, where it was found to be considerably faster and more accurate. Moreover, previous use of Mathieu functions has been limited to constant physical parameters and small degrees of expansions. By contrast, we were able to compute expansions in thousands (and even tens of thousands) of Mathieu functions by making use of the Bessel function expansion of Mathieu–Hankel functions and their asymptotics. This allows quick and robust testing of physical parameters and variations, which may have use in other scattering problems beyond those considered here.

We found that the method coped well with a broad range of frequencies (typically needing more terms for larger k₀, as expected) and smoothly varying porosity/elasticity (with more collocation points and Mathieu functions needed to capture the case of more oscillatory parameters). Our solution representation also directly provides a sine series approximation of the far-field directivity and, unlike standard boundary methods, can easily be evaluated near or on the scatterers. This means that the acoustic near field can be computed efficiently and in a stable manner. These advantages assert that the present method is particularly good for a simple model of turbulence using Lighthill’s analogy.

Examples of diffraction by elastic plates of varying stiffness were presented. We found that a plate with varying stiffness can exhibit an acoustic black hole type behaviour. This has a drastic effect on the near and far fields, both in the scattering and the aeroacoustics setting. Further work is needed to understand how this might be employed as a leading or a trailing edge adaptation to an aerofoil. There is also a potential to use this acoustic black hole effect to move the vibrations away from the trailing edge and into the centre of an aerofoil where they can be baffled.

Finally, we demonstrated that the numerical method can be used on multiple, arbitrary positioned plates. Future work will also look at fast multipole methods and hierarchical solvers for multiple plates and evaluation of the solutions. The method also offers considerable flexibility in the choice of forcing term. In this article, we only considered plane waves and quadrupoles. The new method can easily be extended to boundary conditions different to those in §2 (such as linking different parts of the scatterer or integral constraints) and can be generalized to include boundary conditions on ellipses. While the method is currently restricted to finite plates in two dimensions, it may also be possible to consider similar approaches to other problems (e.g. three dimensions) through separation of variables and different special functions accompanied by spectral methods.

Data accessibility

This work does not contain any experimental data, and all of the results can easily be generated from the equations provided in the article. Numerical code is available at https://github.com/MColbrook/MathieuFunctionCollocation.

Authors' contributions

M.J.C. derived the mathematical model and its solution, developed the numerical method and code, and developed and analysed the examples. A.V.K. developed and analysed the examples in §4 and 5. Both authors contributed to the writing of the manuscript, tested the code and gave final approval for publication.

Competing interests

We declare we have no competing interests.

Funding

This work was supported by EPSRC grant no. EP/L016516/1 (M.J.C.) and Royal Society Dorothy Hodgkin Research Fellowship (A.V.K). The authors thank the Isaac Newton Institute for Mathematical Sciences, Cambridge, for support and hospitality during the programme ‘Bringing pure and applied analysis together via the Wiener–Hopf technique, its generalizations and applications’ where some of the work on this article was undertaken (supported by EPSRC grant no. EP/R014604/1). M.J.C. is also grateful for discussions with Lorna Ayton and Justin Jaworski, and to André Cavalieri and William Wolf for the provision of boundary element code. A.V.K is grateful for useful discussion with Justin Jaworski and Angelis Karlos.