Diffraction of acoustic waves by multiple semi-infinite arrays^a)

Analytical solutions for acoustic/electromagnetic wave scattering problems by different combinations of finite, infinite, and semi-infinite plates/arrays (Lawrie and Abrahams, 2007) are of special interest. These are difficult problems, and we will briefly review some of the work on this subject. Generalising the famous Sommerfeld's half-plane problem (consisting of one semi-infinite plate), Heins (1948a,b) uses the Wiener-Hopf (WH) technique to find the exact solution for the problem where an electric-polarised wave is incident on a pair of parallel semi-infinite plates, symbolising a receiving and transmitting antenna for parts I and II, respectively. This WH formulation has been extended to a matrix version in Abrahams and Wickham (1988, 1990a,b), for the equivalent acoustic problem with a pair of staggered plates, and in Jones (1986) for three semi-infinite planes. Also, by employing appropriate mappings/transformations, it was possible to use the matrix WH techniques for wedges (Daniele and Zich, 2014; Nethercote , 2020a; Shanin, 1998). Additionally, since finding exact solutions is difficult for the matrix WH technique, there have been many developments on approximate and asymptotic factorisation techniques (Kisil, 2018; Kisil and Ayton, 2018; Rogosin and Mishuris, 2016). Alternatively, researchers resort to numerical and asymptotic schemes in their models [see Adams (2008), Craster (2009), Kirby (2008), and Peake and Cooper (2001) for waveguides and ducts examples].

One of the advantages of the WH technique is that a solution for a scattering problem can also be used in aeroacoustics setting to study the interaction of plates with gusts. In particular, Peake (1992) considered an infinite staggered cascade of finite thin blades which were aligned with a uniform subsonic mean flow, and found an iterative solution by an infinite sequence of coupled, semi-infinite WH problems in the high frequency limit. More recently, this work has been extended for non-aligned mean flow (Peake and Kerschen, 1997, 2004), thin aerofoils (Baddoo and Ayton, 2018, 2020), and for high staggering angles (Maierhofer and Peake, 2020, 2022).

In an elastic setting, the scattering and localisation of flexural waves on an elastic Kirchhoff plate is another important problem, especially determining waves that are blocked or trapped (Haslinger , 2016; Haslinger , 2014; Movchan , 2009). Jones (2017) considered a pair of parallel semi-infinite gratings of rigid pins on a Kirchhoff plate. A follow-up article (Haslinger , 2018) studied problems with configurations of parallel semi-infinite gratings. In particular, it featured a problem where four parallel semi-infinite gratings form a waveguide in a herringbone pattern. Both these articles form and use the solution to the discrete WH functional equation but also identify trapped modes from analysing the kernel. The latter of the two articles also assumes that for each side of the waveguide, the two gratings are closely spaced which allowed them to use a dipole approximation.

These diffraction problems can also be considered on a lattice governed by a discrete Helmholtz equation. For instance, waves diffracted in a square lattice by a semi-infinite crack or a semi-infinite rigid constant was solved recently using a discrete WH technique in Sharma (2015a,b), respectively. These two problems are analogous to the classic Sommerfeld's half-plane problem (sound-hard and sound-soft, respectively). This work has been extended to two staggered semi-infinite cracks (Maurya and Sharma, 2019) which was not solved exactly but asymptotically due to the notorious difficulties in factorising a matrix WH kernel.

In this article, we are interested in arbitrary arrays of small Dirichlet cylinders within a continuum. This means that we use the continuous Helmholtz equation subject to boundary conditions imposed on a discrete set of scatterers. The semi-infinite array (the analogue of Sommerfeld's half-plane) problem was solved using the discrete WH technique long ago (Hills and Karp, 1965; Linton and Martin, 2004). However, there has been little work to generalise this to other configurations of arrays as was the case in the continuous and the discrete case outlined above. In our previous work (Nethercote , 2022a), we combined two semi-infinite arrays to form a wedge. Unlike the continuous boundary wedge (Nethercote , 2020a), this did not lead to a matrix WH problem due to the difficulties in finding an appropriate mapping. Instead, we considered two coupled systems of equations which were solved using the discrete WH technique, followed by an effective numerical iterative procedure. In this article, we will study problems involving any number of independent semi-infinite arrays comprised of equidistant scatterers. This is very general since the position and orientation of the arrays are arbitrary, which allows us to model many types of interesting problems. As in Nethercote (2022a), we will use the WH technique for each of the arrays and then couple them together. Yet this time, the coupling is encoded directly in the matrix inversion which means that there is no need for an iterative scheme.

For any arbitrary configuration of scatterers, one can use numerical techniques to determine the scattering behaviour. Some examples of these techniques include finite element methods, a T-matrix reduced order model (Ganesh and Hawkins, 2018; Hawkins, 2023) and a least square collocation approach that was used in Chapman (2015) and Hewett and Hewitt (2016) to study the electrostatic and electromagnetic shielding by Faraday cages. While these methods are very efficient at modelling the interactions between individual scatterers, they do not work very well at modelling the infinite nature of periodic arrays.

The structure of the paper is as follows. We start by setting up and solving the WH problem in Sec. II A, which results in a matrix equation that is inverted in Sec. II B to find the scattering coefficients. In Sec. II C, we proceed by looking into the special case with two semi-infinite arrays and link with the point scatterer wedge from Nethercote (2022a). We also analyse the determinants and condition numbers of the matrices involved in Sec. II D as well as discuss the use of fast multipole methods for efficient computations of their components in Sec. II E. Finally, we showcase several different test cases in Sec. II F and compare the WH solution with other numerical techniques in Sec. II G and highlight their strengths and weaknesses.

Viewed as a three-dimensional problem, the scatterers are all cylinders of infinite height, have a small radius, and satisfy homogeneous Dirichlet boundary conditions. This problem can naturally be reduced to two dimensions for non-skew incidence and this is what we will be considering here. Throughout this article, we will exploit the methodology described in Nethercote (2022a), since the point scatterer wedge (or line scatterer wedge when viewed in three dimensions) can be considered as a particular case of what is presented here.

It is fairly simple to match these results with the specific case of the point scatterer wedge studied in Nethercote (2022a). That wedge configuration is produced from these parameter choices: $a 1 = a 2 = a , s 1 = s 2 = R 0 ( 2 ) = s , R 0 ( 1 ) = 0 , θ 0 ( 2 ) = α 2 = − α 1 = − α$ which gives the distance function $Λ ( 1 , 2 ) ( m , n ) = s ( m 2 + ( n + 1 ) 2 − 2 m ( n + 1 ) cos ( 2 α ) ) 1 / 2$⁠. This leads to the two WH kernels [ $K 1 ( z )$ and $K 2 ( z )$] as well as the resulting coefficients (⁠ $λ 1 , n$ and $λ 2 , n$⁠) becoming identical. There are two key differences left. One is that the second array here has one extra scatterer after truncation, i.e., the vector $A ( 2 )$ has one extra entry. However, this extra entry can be neglected and consequently, one must ignore the last column in $M ( 1 , 2 )$ and the last row in $M ( 2 , 1 )$ as well. The other difference is that we have revoked the need for the iterative scheme and have instead inverted the matrix equation (32) directly. Note that the inverted matrix can be expanded, $( I − M ( 1 , 2 ) M ( 2 , 1 ) ) − 1 = I + M ( 1 , 2 ) M ( 2 , 1 ) + ( M ( 1 , 2 ) M ( 2 , 1 ) ) 2 + ⋯$⁠, which can then be used to recover every iteration in the iterative scheme [although this requires the spectral radius $ρ ( M ( 1 , 2 ) M ( 2 , 1 ) ) < 1$ to converge].

The behaviour of the full matrix is naturally very different because the condition number for the matrix in Eq. (30) is quite moderate (approximately $10 0 − 2$⁠) and the determinants are non-zero for all cases that we tested for this article. If the system had a zero determinant, it could either be due to the modelling assumption being insufficient (Nethercote , 2022b) or due to a specific physical phenomenon. Since a zero determinant implies a non-unique solution to the matrix equation and considering that this mostly relies on the chosen geometry of the problem, one could conjecture that this could imply the presence of homogeneous (Rayleigh-Bloch) waves. However, since it was proven that semi-infinite arrays with Dirichlet boundary conditions cannot support these waves (Bonnet-Ben Dhia , 2016; Bonnet-Ben Dhia and Starling, 1994), we are not expecting a zero determinant.

To calculate the matrices $M ¯ ( j , ℓ )$⁠, we can use the fast multipole methods (FMM) library which is accessible from Greengard and Gimbutas (2012) [see also Beatson and Greengard (1997) for details on the algorithm]. This method is very accurate at computing sums of the form (35). For the matrices $M ( j , ℓ )$⁠, the algorithm is able to reduce the computational cost with respect to the truncation (N) from $O ( N 3 )$ to as low as $O ( N 2 ln ( N ) )$⁠. However, it is only able to do this if the values of ks_j are sufficiently small for all j. To demonstrate this, Fig. 4 plots the computation times to calculate the matrix in Eq. (30) for the point scatterer wedge given in Fig. 5(a) (where $s 1 = s 2 = s$⁠). This figure shows the difference between both methods with respect to the truncation N and ks. On the left side, we see that for smaller ks the computational order has been reduced. On the right side, we see that FMM becomes slower than using direct methods when ks is approximately larger than π. Although there are developments for larger wavenumbers (Crutchfield , 2006), this has not been implemented in the library we used for the current work. As a result, we will only use FMM to calculate $M ¯ ( j , ℓ )$ if the values of ks_j are sufficiently small for all j. If not, then we will calculate it directly.

In this section, we consider and showcase several different examples of test cases in Fig. 5 by plotting the real part of the total wave field. We look at five different configurations of semi-infinite arrays, where the array parameters are given in Table I.

The first case [Fig. 5(a)] is an exact recreation of the same point scatterer wedge that was considered in Fig. 6 of Nethercote (2022a). This time, however, we are able to create gaps in the wedge interface [Fig. 5(b)] as well as add additional scatterers to one or both of the arrays [Fig. 5(c)], provided that we do not have overlapping scatterers. With this new generalised formulation, we can also consider cases with additional arrays. For example, we can have twelve outwardly pointing arrays where the ends are positioned to create a cage [Fig. 5(d)]. This cage is able to shield the middle region from an incident wave with a low wavenumber. In particular, the sound pressure level inside the cage [given by the formula $20 log 10 ( root mean square ( Φ ) )$] is approximately –26.36 compared to 0 when there are no scatterers at all. This configuration is analogous to electrostatic and electromagnetic shielding problems using Faraday cages [see Chapman (2015) and Hewett and Hewitt (2016)].

Although the WH technique is not typically used outside of semi-infinite problems, it is important to model the wave scattering by other types of arrays including circular arrays (Martin, 2014), infinite arrays with defects (Thompson and Linton, 2008), and long finite arrays (Thompson , 2008). Another case of special interest is determining the bandgap structure of doubly periodic lattices (Botten , 2001). For example, McIver (2007) and Krynkin and McIver (2009) study an infinite doubly periodic lattice of small scatterers with Dirichlet boundary conditions. With that in mind, the final case that we consider is a series of stacked infinite arrays to create a finitely thick doubly periodic lattice. The reason why we chose to look into this case is because we wanted to know if this configuration has similar properties to the fully periodic lattice as told by the bandgap diagrams of (Krynkin and McIver, 2009). Specifically, we wanted to see the behaviour of an incident wave that cannot penetrate the lattice (i.e., in a stop band) and the Bloch waves resulting from one that can (i.e., in a passband). Choosing a wavenumber within a stop band [Fig. 5(e)] causes the incident wave to be almost fully reflected from the lattice. The alternative choice of a wavenumber within a passband [Fig. 5(f)] does cause some reflection but most of the energy goes into the lattice and forms a Bloch wave inside before becoming a transmission out of the other side.

Here, we seek to find a means of comparison with numerical methods that is both as fair and comprehensive as possible. We have three possible methods to compare with: finite element software comsol, a T-matrix solver (TMAT) by Hawkins (2023) and a least square collocation (LSC) method which was used in Chapman (2015) and Hewett and Hewitt (2016) for solving Laplace's equation and Helmholtz's equation, respectively.

In previous articles (Nethercote , 2020b,a), we have extensively used comsol for comparison. However, there are limitations to the comparison since comsol is not able to find a solution with thousands of scatterers to fully compare with our method. Here, we will compare with results computed using TMAT or LSC where we are able to have the same number of scatterers and also the computation is performed on the same computer. We are also able to compute the scattering coefficients with these methods which is not possible with comsol.

The T-matrix software package provides an object-oriented implementation of an efficient reduced order model framework for modelling two- and three-dimensional wave scattering simulations. For the scattered field, TMAT uses the multipole expansion which is truncated depending on the scatterer size relative to the wavelength. TMAT also creates a Bessel function expansion of the incident field about every single scatterer which is truncated to the same number of terms as the multipole expansion. Following the construction of the T-matrix for every scatterer, a system of matrix equations is formed and solved for the scattering coefficients. However, the current version of this software restricts the truncation such that dipole coefficients must be included as well as monopole coefficients.

The idea of the least squares collocation method is to create and solve an overdetermined matrix system to find the scattering coefficients of a truncated multipole expansion. Here, the known data is a collection of boundary data from a number of collocation points for every scatterer. We would need more collocation points per scatterer than the number of multipole terms to guarantee that the matrix system is overdetermined. The TMAT and LSC methods are both able to consider a large number of scatterers and despite the different methods, they are comparable in their results such that we only need to compare the WH method with one of them. We choose the LSC method because it will be a fairer comparison since we can restrict that one to monopole coefficients only.

With this in mind, Fig. 6 is a collection of plots of the absolute value difference between the two sets of scattering coefficients where the truncation is chosen to be 1000. In these plots, the index n of the coefficients is ordered such that $A − 1001 , … , A − 1 , A 0 , … , A 1000$ corresponds to $A 1000 ( 2 ) , … , A 0 ( 2 ) , A 0 ( 1 ) , … , A 1000 ( 1 )$⁠, respectively. On the top row of Fig. 6, we look at the infinite array problem and have an incident wave with wavenumber $k = 5 π$ but two different incident angles; $θ I = π / 4$ and $θ I = π / 12$ for the left and right side, respectively. Note that equivalent plots of the relative error share the same shape and are on similar scale as Fig. 6. On the bottom row, we look at the point scatterer wedge given in Table I and use the same incident waves as in Fig. 6 of Nethercote (2022a) (i.e., $k = 5 π$ and $θ I = 0$ for the left side and $k = 15 π$ and $θ I = π / 2$ for the right side).

The top row of Fig. 6 is especially interesting because the comparison with the exact solution allows us to decompose the error between the WH and LSC methods and assess their strengths and weaknesses. One conclusion of this decomposition can be seen in the middle of the plots (⁠ $n ≈ 0$⁠) where the WH seems to be at its weakest and the LSC method at its strongest. This is due to how each method sees the problem as the WH method considers each array separately and adds the interaction between arrays when solved, whereas the LSC method considers each scatterer separately and solves for the individual interactions. The other conclusion of this decomposition can be seen at the truncated ends of the arrays (⁠ $n ≈ ± 1000$⁠) where the WH is at its strongest and the LSC method at its weakest. This is again due to how the methods see the problem as separate semi-infinite arrays or individual scatterers. It is likely that we will come to similar conclusions for most (if not all) of the potential configurations of this setup. It is also likely that the weak region for the WH method will improve when the semi-infinite arrays are well separated. It is important to note that the overall error does converge to zero as the truncation increases and shape of these error graphs scales with the truncation as well.

To summarise, we have generalised the WH method used for diffraction by a wedge of point scatterers (Nethercote , 2022a), by considering any number of semi-infinite arrays with an arbitrary set of parameters. The method remains essentially unchanged with the additional benefit of having removed the need for an iterative scheme. The matlab scripts we created from this solution are quite versatile and we are able to consider a very wide range of cases, some of which are illustrated in Fig. 5.

We have also compared the WH method with some numerical approaches such as the LSC method, and Fig. 6 has shown that the two methods have some good agreement and highlighted the strengths and weaknesses between them. We found that the LSC method is better at modelling the interactions between the scatterers at the ends of the arrays and the WH method is better at modelling the infiniteness of the arrays. Knowing this, one could propose to use a hybrid of the two methods to get accurate coefficients for all scatterers (in other words, use LSC to get the coefficients $A m ( j )$ for small m and WH for large m). In theory, this hybrid would have the strengths of both methods but neither of the weaknesses. While finding the optimal m where we should switch methods would be simple for infinite array cases, more general configurations will be more difficult. This is because we will not have an exact solution to decompose the error quantity $| WHT − LSC |$ and this optimal m will not be unique.

Although we have not explicitly discussed resonance in the current article [see Nethercote (2022b) for an overview], it is nonetheless of special interest. By exclusively using the methods discussed in this article, we are capable of numerically evaluating cases where inward resonance is occurring, $( k s j / 2 π ) ( 1 + cos ( α j − θ I ) ) ∈ ℤ$⁠. However, we cannot numerically evaluate outward resonance (and by extension double resonance) cases, $( k s j / 2 π ) ( 1 − cos ( α j − θ I ) ) ∈ ℤ$⁠, because the associated scattering coefficients will tend to zero but still lead to a non-trivial wave field. Finding a general procedure which can find and extract outward resonant waves is a topic for future work.

Another interesting avenue to pursue is to find a way to use the solution for the scattering coefficients to identify special features. Examples of this could include the parameters of Bloch waves in lattices [see the multiple infinite array case given by Fig. 5(f)], trapped modes in waveguides or the constructive/destructive interference caused by waveguides or the acoustic Faraday cage.

This research was supported by EPSRC grant EP/W018381/1. A.V.K. is supported by a Royal Society Dorothy Hodgkin Research Fellowship and a Dame Kathleen Ollerenshaw Fellowship. A.V.K. also acknowledges the support from EU through the H2020-MSCA-RISE2020 project EffectFact, Grant agreement ID No. 101008140. M.A.N. was supported by a David Crighton fellowship at Cambridge University and thanks Nigel Peake for many insightful discussions which greatly refined this article. The authors would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme Mathematical theory and applications of multiple wave scattering when work on this paper was undertaken. This programme was supported by EPSRC Grant No. EP/R014604/1. The authors also give thanks to Stuart Hawkins and David Hewett for providing the source code for the TMAT and LSC methods, respectively.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.