Analytical methods for perfect wedge diffraction: a review
Matthew A. Nethercote, Raphael C. Assier, I. David Abrahams

TL;DR
This review surveys various analytical methods for solving the classical problem of wave diffraction by a wedge, comparing their effectiveness and introducing some novel approaches.
Contribution
It compiles and compares multiple analytical techniques for wedge diffraction, including less-known methods, and provides a critical assessment of their accuracy and practicality.
Findings
Exact solutions are obtained using multiple methods.
GTD approximations are compared with exact solutions.
New methods like the embedding approach are introduced.
Abstract
The subject of diffraction of waves by sharp boundaries has been studied intensively for well over a century, initiated by groundbreaking mathematicians and physicists including Sommerfeld, Macdonald and Poincar\'e. The significance of such canonical diffraction models, and their analytical solutions, was recognised much more broadly thanks to Keller, who introduced a geometrical theory of diffraction (GTD) in the middle of the last century, and other important mathematicians such as Fock and Babich. This has led to a very wide variety of approaches to be developed in order to tackle such two and three dimensional diffraction problems, with the purpose of obtaining elegant and compact analytic solutions capable of easy numerical evaluation. The purpose of this review article is to showcase the disparate mathematical techniques that have been proposed. For ease of exposition,…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Analytical methods for perfect wedge diffraction: a review
Matthew A. Nethercote*∗, Raphael C. Assier∗* and I. David Abrahams*†*
∗ School of Mathematics, University of Manchester, Oxford Road, Manchester, M13 9PL, UK
† Isaac Newton Institute, University of Cambridge, 20 Clarkson Road, Cambridge CB3 0EH, UK
Abstract
The subject of diffraction of waves by sharp boundaries has been studied intensively for well over a century, initiated by groundbreaking mathematicians and physicists including Sommerfeld, Macdonald and Poincaré. The significance of such canonical diffraction models, and their analytical solutions, was recognised much more broadly thanks to Keller, who introduced a geometrical theory of diffraction (GTD) in the middle of the last century, and other important mathematicians such as Fock and Babich. This has led to a very wide variety of approaches to be developed in order to tackle such two and three dimensional diffraction problems, with the purpose of obtaining elegant and compact analytic solutions capable of easy numerical evaluation.
The purpose of this review article is to showcase the disparate mathematical techniques that have been proposed. For ease of exposition, mathematical brevity, and for the broadest interest to the reader, all approaches are aimed at one canonical model, namely diffraction of a monochromatic scalar plane wave by a two-dimensional wedge with perfect Dirichlet or Neumann boundaries. The first three approaches offered are those most commonly used today in diffraction theory, although not necessarily in the context of wedge diffraction. These are the Sommerfeld-Malyuzhinets method, the Wiener-Hopf technique, and the Kontorovich-Lebedev transform approach. Then follows three less well-known and somewhat novel methods, which would be of interest even to specialists in the field, namely the embedding method, a random walk approach, and the technique of functionally-invariant solutions.
Having offered the exact solution of this problem in a variety of forms, a numerical comparison between the exact solution and several powerful approximations such as GTD is performed and critically assessed.
1 Introduction and formulation
At the close of the century, wedge diffraction became a core problem in mathematical physics when renowned mathematicians Poincaré and Sommerfeld studied the diffraction of wave fields in angular domains (Poincaré, 1892, 1897; Sommerfeld, 1896, 1901). Sommerfeld made the first breakthrough when he solved his famous half-plane problem (Sommerfeld, 1896), during which he introduced the contour integral representation that we now know as the Sommerfeld integral. This work has now been translated to English in Sommerfeld (2003), with additional insightful comments. Sommerfeld would later be the first to solve problems of wedge diffraction (Sommerfeld, 1901) where the wedge has an interior angle equal to ().
For wedges with arbitrary interior angles, the solution was first obtained by Macdonald (1902). He did this by considering a line source incident wave and used separation of variables to get a series solution. The solution was rewritten in Sommerfeld integral form and he then provided the solution for an incident plane wave. We discuss this line source approach and provide an alternative way to obtain the plane wave solution in A.
In the 1950s, Malyuzhinets released a series of papers that culminated in the solution to the problem with impedance boundary conditions, (Malyuzhinets, 1955a, b, 1958b, 1958c, 1958a). This result created the first method that we discuss here in Section 2, the Sommerfeld-Malyuzhinets technique (S-M). Other authors who solved the impedance wedge problem independently were Senior (1959) and Williams (1959) but for more details on Malyuzhinets’ method, see the review paper (Osipov and Norris, 1999) or the books (Budaev, 1995; Babich et al., 2007; Lyalinov and Zhu, 2013).
One of the most popular methods in diffraction theory is the Wiener-Hopf (W-H) technique, invented by Wiener and Hopf (1931) as a means to solve a special type of integral equation. It was soon discovered to be a useful method for diffraction problems and has appeared in a number of classic articles such as (Copson, 1946) and (Jones, 1952). Since then, applications of the technique have appeared in a wide array of research areas including diffraction, waveguides and flow problems.
The well-known textbook (Noble, 1958) provides an excellent tutorial for various aspects and extensions of the W-H technique. In 2007, the Journal of Engineering Mathematics published a W-H special issue led by a historical overview (Lawrie and Abrahams, 2007) along with a collection of articles applying the W-H technique to various problems. For wedge problems, the technique was thought to be ineffective due to the two boundaries not being parallel, however (see e.g. (Shanin, 1996; Daniele, 2003b)) this can be overcome as discussed later in Section 3.
Another key method is based on the Kontorovich-Lebedev (K-L) transform. First introduced by Kontorovich and Lebedev (1939), this transform is an effective tool when dealing with a radial coordinate. This makes it useful for wedge diffraction problems as evidenced by Abrahams (1986, 1987), since obtaining the general solution in that way is a very natural process. We shall discuss this further in Section 4, but for more details on the transform and its applications, see (Lebedev, 1965), (Jones, 1964, 1980, 1986) and (Felsen and Marcuvitz, 1994).
In Section 5, we will focus on the asymptotic technique created by the classic paper (Keller, 1962), called the Geometrical Theory of Diffraction (GTD), see also the book by Borovikov and Kinber (1994). We will also follow the uniform GTD extension detailed in literature such as (Kouyoumjian and Pathak, 1974; James, 1986; Mcnamara et al., 1990) and (Babich et al., 2007).
Even though we will not use it in this review, for completeness, it is important to mention an alternative asymptotic technique applied to diffraction problems that is the Physical Theory of Diffraction (PTD) (Ufimtsev, 1971). This development was made possible in part thanks to Macdonald’s work on Kirchhoff’s approximation (see e.g. Ufimtsev (2014)). A useful paper that compares the GTD and PTD asymptotic techniques as well as the exact solution in series and integral form is (Hacivelioglu et al., 2011). Similar methods, describing creeping waves in diffraction by smooth obstacle, have also been developed in Fock (1965) for example.
Section 6 contains a number of alternative methods that are effective but less well-known for wedge diffraction. The first of these alternative methods is based on the very powerful concept of embedding formula. This reasonably recent approach consists in expressing the diffraction coefficient (which depends on both the incident and observer angles) of the diffracted field resulting from an incident plane wave in terms of the directivities (depending on one angle only) of simpler problems. These simpler problems are directly related to edge Green’s functions. These are Green’s functions for which the source is sent towards the geometric singularities of the obstacle. This method was primarily used for planar cracks and slits, and parallel combinations of these (see e.g. (Williams, 1982; Gautesen, 1983; Martin and Wickham, 1983; Biggs, 2001, 2002)). In (Craster and Shanin, 2005) it was shown that the method can be successfully adapted to wedges, as we will discuss later.
The second of these alternative methods is the so-called random walk approach. It is based on the known link between deterministic PDEs and stochastic differential equations (SDEs) given by the Feynman-Kac theorem. It allows to express the solution of a diffraction problem as the mean of a set of solutions to given SDEs with carefully chosen initial and final conditions. The method was developed through a series of papers by Budaev and Bogy (2001, 2002a, 2002b, 2003), the latter being dedicated to wedge diffraction.
The last of these is the method of functionally-invariant solutions also known as the Sobolev-Smirnov method which has been used for a number of plane wave diffraction problems from half-planes (Sobolev, 1935; Smirnov, 1964) to wedges (Filippov, 1964; Komech et al., 2015; Babich, 2015). A very similar method that develops Busemann’s “conical flow method” (Busemann, 1947) was also considered in Keller and Blank (1951) and Miles (1952) for example.
In this review (apart from Macdonald’s approach discussed in A), we will focus primarily on plane wave incidence rather than line sources. It has to be noted, however, that a broad range of work (Bromwich, 1915; Oberhettinger, 1954; Rawlins, 1987, 1989) has also been carried out for both acoustic and electromagnetic sources. For other reviews of some of the methods used for various types of incident waves (plane, cylindrical, spherical, dipole and pulse), see Oberhettinger (1958) and Bowman et al. (1987).
The elastic wedge problem has equally received a lot of attention. Knopoff (1969) wrote an interesting review of possible approaches to tackle this (still unsolved) problem, it includes attempts using the method of images, the W-H technique, the K-L transform and the conical flow method. More recent approaches by Croisille and Lebeau (1999) or Budaev and Bogy (1995, 1996, 1998) are also worth mentioning.
In Section 3.1, in the spirit of Wegert (2012), we will visualise complex functions using phase portraits in order to show domains of analyticity, locations of singularities and orientations of branch cuts. These portraits assign the argument of a complex function to a HSV colour model. For example, Figure 1 (left) shows the phase portrait of which we use as a colour reference.
Wedge diffraction has a number of physical applications such as the scattering of acoustic pressure fields or electromagnetic fields by sharp structures, scattering of the Sun’s radiation by cloud ice crystals and seismology. We shall study this from a more mathematical perspective. From here onward, we will assume that the problem is time-harmonic with time factor , which is therefore suppressed, and we will consider solutions to the homogeneous Helmholtz equation,
[TABLE]
inside a wedge-shaped region described in polar coordinates by (see Figure 1 (right)). The complementary region, is considered to be the wedge scatterer. Defining , the interior angle of the wedge scatterer is . Throughout this paper, we will consider two cases of homogeneous boundary conditions (BCs), Dirichlet and Neumann, on both faces of the wedge,
[TABLE]
For acoustics, (1.2) and (1.3) are called sound-soft and sound-hard BCs respectively. For electromagnetic scattering by an electric (resp. magnetic) polarized plane wave, the solution to the perfect electric conducting (PEC) problem can be expressed in terms of a potential satisfying (1.2) (resp. (1.3)).
We define the incident plane wave as, with wavenumber and incident angle . Due to the symmetry of the problem, we restrict the incident angle to . Figure 1 (right) illustrates the geometry of the problem.
An initial approximation is found using classic Geometrical Optics (GO). The GO part of the solution consists of the incident wave and any reflected waves produced. The rest of the solution is considered to be the diffracted field, which satisfies a two-dimensional radiation condition (see (Schot, 1992) for a good review on this), written in integral form:
[TABLE]
Lastly, there will be an edge (or Meixner) condition as becomes small.
[TABLE]
and for the case of Dirichlet BCs. Typically, the edge conditions can be derived using the Frobenius method (Bender and Orszag, 1999) while ensuring that the energy remains finite in any neighbourhood of the wedge edge.
2 The Sommerfeld-Malyuzhinets technique
The first method to be reviewed is the Sommerfeld-Malyuzhinets (S-M) technique. Here we will show briefly how to get the solution to the perfect wedge problem (1.1)-(1.5). For a more thorough explanation, consult Sections 1-4 of (Babich et al., 2007). This technique is based on the general solution of diffraction problems in angular domains being represented as the Sommerfeld integral
[TABLE]
where are contours defined in Figure 2 (left) and is an unknown function to be determined. This representation ensures that the Helmholtz equation is automatically satisfied111This is proven using integration by parts and noting that satisfies the Helmholtz equation with polar coordinates .. The form of the spectral part, , is necessary for the radiation conditions to be satisfied. This can be proven by using the method of steepest descent to approximate the integral as (see section 3.7 of Babich et al. (2007)). The function , referred to as the spectral function, is assumed to be meromorphic in the domain
[TABLE]
for some and analytic in the same domain with , see Figure 2 (right). The poles of will be seen to correspond to the geometrical optics part of the wave field. The structure of (2.1) means that is defined up to an additive constant. Using the edge conditions, we can assume that has the following behaviour as ,
[TABLE]
where for the Dirichlet case.222We assume (2.3) for convenience later. Say we assumed as where instead. Then a later step will require us to redefine using the additive constant property such that implying (2.3). It is also important to note that is an entire -periodic function of , which decays rapidly as only in the set of half-strips,
[TABLE]
displayed in Figure 2 (left), and grows rapidly in the complementary set. The Sommerfeld contours are defined so that the integrand is analytic and decays as along these contours.333Note that the Sommerfeld contours are contained in the domain (2) with .
A crucial part to the S-M technique is Malyuzhinets’ Theorem or the Sommerfeld Nullification Theorem. This is an important theorem because it allows to obtain a functional equation satisfied by the spectral function. The theorem and its proof are presented by Malyuzhinets in (Malyuzhinets, 1958b) and more recently in Section 3.4 of (Babich et al., 2007).
Theorem 1** **(Malyuzhinets’ Theorem or Sommerfeld Nullification Theorem)
Let the function be analytic and single-valued inside the half-strip,
[TABLE]
If for some constant , the function has the following behaviour as in this half-strip,
[TABLE]
and for any ,
[TABLE]
then,
[TABLE]
where is the integer part of and the coefficients are constants.
With Malyuzhinets’ theorem, we have all the tools required to determine the spectral function . Recall that it was assumed to be meromorphic in the domain (2), moreover it has only one simple pole with unit residue at within the strip correponding to the incident wave .
2.1 Dirichlet boundary condition
Applying the Dirichlet BCs (1.2) to the general solution (2.1) implies that
[TABLE]
Due to (2.3), we can apply Malyuzhinets’ theorem to (2.10) to produce a pair of functional equations for ,
[TABLE]
These equations imply that is symmetric about and as a consequence is periodic. Because the pole at produces the incident wave, its residue should be 1, i. e.
[TABLE]
Because this pole is the only singularity in the strip , then by the determined symmetry, we also have poles at and with residue . These two poles correspond to the top and bottom reflected waves respectively. The periodicity implies that each of the poles are repeated every with the same residue. We can therefore express as a sum of poles,
[TABLE]
Using the definition and pole expansion of ,
[TABLE]
we rewrite as follows,
[TABLE]
It is easy to check that (2.13) satisfies the functional equations (2.11), satisfies the estimate as , and has a single pole with unit residue in the strip . This means that the solution to the exterior wedge problem with Dirichlet BCs is,
[TABLE]
2.2 Neumann boundary condition
Solving for the Neumann case is done in a very similar way. Applying the Neumann BCs (1.3) to the general solution (2.1) implies that
[TABLE]
We integrate by parts and apply the Malyuzhinets’ theorem to (2.15) to obtain the following pair of functional equations for ,
[TABLE]
Applying (2.3) determines that . Then the functional equations become,
[TABLE]
These equations imply that is antisymmetric about . This symmetry also shows that is a periodic function. Because is the only pole of in the strip and has unit residue, then due to the antisymmetry and periodicity of , there are more poles located at and for , all with unit residue. We can therefore express as a sum of poles,
[TABLE]
It is easy to check that (2.18) satisfies the functional equations (2.17), satisfies the estimate (2.3) as , and has one simple pole with unit residue in the strip . Finally, the solution to the exterior wedge problem with Neumann BCs is,
[TABLE]
Critical analysis
The main advantage of the S-M technique is its ease of implementation. Once in the right form, a natural deformation to the steepest descent contour can transform it into a very simple and fast converging integral (see Section 5). This aspect clearly makes it the gold standard representation of the general solution to the perfect wedge problem. Moreover, because of the form of the integrand, one can somehow think of the solution as a weighted superposition of plane waves.
An other advantage of such method is its flexibility. It is indeed possible to use it for more complicated problems such as the impedance wedge (Malyuzhinets, 1958a) or with various types of incident waves (Bowman et al., 1987). It has also been used in the quarter-plane problem to infer some results about the far-field structure (Lyalinov, 2013).
A disadvantage of the method is that it is not particularly constructive. In fact, it often starts with a form of the solution being posed, which, may somehow appear as some kind of “black magic”. The B addresses this issue by showing that the two contour representation of the solution comes naturally from Green’s identity. Moreover, once the integral form of the solution is written down, it is not straightforward to prove that the radiation condition is satisfied.
3 The Wiener-Hopf technique
The second method to be reviewed is the Wiener-Hopf (W-H) technique. Before authors such as Shanin and Daniele used the W-H technique for wedge problems, it was mostly used for waveguide problems or more complicated half-plane problems. In the two papers (Shanin, 1996, 1998), Shanin looks at solving wedge problems with inhomogeneous impedance BCs via the W-H technique. In a series of papers and internal reports (Daniele, 2000, 2001, 2003a, 2003b), Daniele develops several aspects of his method for impenetrable wedge problems.
In this section, we will combine elements from both (Shanin, 1996) and (Daniele, 2003b) to rederive (2.14) and (2.19) using the W-H technique. The idea is to Laplace transform and its derivative on the two wedge faces, , and the line of symmetry . These transforms are used to produce the W-H equations. After the BCs are considered, a mapping to a new complex plane is introduced so that the W-H technique can be applied to produce a solution that will match (2.14) and (2.19). We will start by studying this mapping.
3.1 A useful mapping
As we will see later, when the W-H equations are derived, they cannot be easily factorised using standard methods. To counter this issue, we will need to map these equations onto a new complex plane so that they can be reduced to classical W-H equations like those in (Noble, 1958). In order to do that, Shanin and Daniele use slightly different mappings,
[TABLE]
where and are the old and new complex variables. Though the two mappings are conceptually equivalent, we will study Daniele’s mapping in what follows. Note though that in his work, Daniele has assumed that has a small imaginary part in order to have a strip of analyticity for the W-H equations. However this is not strictly necessary, and here we will consider , essentially reducing the W-H problem to a Riemann-Hilbert problem (see e.g. (Kisil, 2015)). This means that does not need to appear explicitly in the mapping, and we can simply use
[TABLE]
with the corresponding inverse,
[TABLE]
We also consider the intermediate mapping and corresponding inverse,
[TABLE]
The mapping (3.1) has a single branch cut along the real line segment where the local argument of the chosen branch is . This is done by choosing the branch of the inverse cosine such that , which is standard for programs such as Mathematica and MATLAB. Note that the intermediate mapping limits to belong to the strip .
One of the most important features of the mapping (3.1) is that the upper half plane (UHP) is mapped to a subset of the UHP as shown in Figures 3a and 3b. This implies in particular that if a function is analytic in the UHP, then the function is analytic in the UHP. Another noteworthy property is that if a function has no branch point at (or ), then the function will be symmetric with respect to (or ) in the -plane.
The mapping (3.1) is designed specifically for manipulation of the following functions,
[TABLE]
where the branch for the square root is chosen such that . For context, is used to identify domains of analyticity whereas, and are kernel functions that need to be factorised. Noting that and have the relation , we map to the and planes.
[TABLE]
In the plane, has branch cuts on the segments and (see Figure 4b which is a phase plot of on the -plane).
Similarly we can map in the and planes.
[TABLE]
In the plane, has a single branch cut on the segment and is related to by the identity , as shown by comparing the phase portraits in Figures 4a and 4c. The consequence is that the lower half plane (LHP) is mapped to a subset of the region as can be seen from comparing Figures 3c and 3d. This also implies, in particular, that if a function is analytic in the UHP, then the function is analytic in the LHP.
Lastly, we study the effect of the mapping on ,
[TABLE]
showing that and are closely related in the -plane, in the sense that and have the same branch cuts but their phase plots are the symmetric images of each other about , as illustrated in Figures 4b and 4d. This relationship means that factorising in the plane is analogous to factorising .
Before performing such factorisation, let us define the domains on which this factorisation will be done:
[TABLE]
Note that and . Factorising , requires to write
[TABLE]
where is analytic in and is analytic in . We expect that will have a branch cut starting at and a branch cut starting at . Hence and will be symmetric about the points and respectively. We can realise this factorisation by using the fact that both the leading order behaviour of as and the jump across the cut are consistent with the function . This means that we can define and by,
[TABLE]
Clearly has a branch cut on the segment and is analytic at . While retains the branch cut on the segment , and dividing by has made become a removable singularity and cancelled out the cut discontinuity. Hence we can assign the limiting value to and make analytic at that point. Recalling the definition , we also map and to the -plane as follows:
[TABLE]
As anticipated, due to the absence of the branch point at , is symmetric about . Similarly, due to the absence of the branch point at , is symmetric about . Figure 5a illustrates this factorisation by the relationship (3.10).
Now we factorise using the established relation (see (3.7)) to find that and . It can be shown that has a branch cut on the segment , is analytic at and is symmetric about , while has a branch cut on the segment , can be made analytic at and is symmetric about .
3.2 Derivation of the Wiener-Hopf equations
Daniele derives the W-H equations by rewriting the Helmholtz equation in terms of an oblique Cartesian coordinate system and uses a Laplace transformation in each of the new coordinates (Daniele, 2003b). However this process can be time-consuming and we will show a different method here.
(Shanin, 1996) tackles an interior wedge problem with inhomogeneous impedance BCs. In this geometry, only one W-H equation is derived, which is obtained via Green’s second identity. However we will need to split the exterior wedge region into two halves to obtain two W-H equations. Take Green’s second identity for functions , twice continuously differentiable on domain with boundary ,
[TABLE]
where is the normal derivative. Here is the unknown solution and we choose a suitable test function for that satisfies the Helmholtz equation (1.1). Then the left hand side of (3.13) is automatically zero. We do this for two wedge regions and which require different test functions. The right hand side of (3.13) has two parts, the wedge boundary at and an imaginary boundary at . For the upper region, we choose the test function , leading (3.13) to become
[TABLE]
For the lower region, , we choose the slightly modified test function leading (3.13) to become
[TABLE]
Define the Laplace transform with the following inverse,
[TABLE]
where is analytic in the half-plane , then we define the transforms of and ,
[TABLE]
These transforms are applied to both (3.14) and (3.15) which produces the W-H equations,
[TABLE]
Adding and subtracting these two equations leads to,
[TABLE]
These are the so-called generalised W-H equations. In this system, the functions and are analytic in a region containing the upper half plane regardless of the value of . We solve the system (3.19) for and using the boundary data on the right hand side. Noting that is equal to with the opposite orientation, the inverse transform of is,
[TABLE]
which is clearly very similar to the Sommerfeld integral. Applying Malyuzhinets’ theorem, we find that,
[TABLE]
We can derive a second formula by comparing the inverse transform of with the following,
[TABLE]
which is obtained by differentiating (2.1) with respect to and integration by parts. Applying Malyuzhinets’ theorem, we find that,
[TABLE]
Adding (3.21) to (3.23) and setting implies that for all ,
[TABLE]
establishing a link between the spectral function and the W-H unknowns. We will now apply the BCs to solve the W-H system (3.19).
3.3 Dirichlet boundary condition
The transformed Dirichlet BCs are which simplify (3.19) to,
[TABLE]
In this form the W-H technique cannot be applied, so (3.3) and (3.6) (discussed in Section 3.1) are used here,
[TABLE]
In the -plane, and are analytic in , except for some potential poles on the real line segment . Similarly, are analytic in , except for some potential poles on . We have already factorised in (3.11), leading (3.27) to become,
[TABLE]
For both equations, (3.26) and (3.28), the left sides are meromorphic in and the right sides are meromorphic in , however due to the potential poles, these equations cannot be used to create an entire function. To counteract this we remove the poles on the right side using the knowledge of the GO component of the solution. Assuming that , the GO components of are,
[TABLE]
where is the Heaviside function. The two poles correspond to , i.e. to . However there is no in the chosen branch of inverse cosine that satisfies . This means that the only pole that needs to be removed is that of at , corresponding to . The residue at this pole is,
[TABLE]
With this residue, we remove the pole from equations (3.26) and (3.28),
[TABLE]
In both equations (3.30) and (3.31), the left sides are now analytic in and the right sides are analytic in . In order to apply Liouville’s theorem, we must determine the behaviour of each part in equations (3.30) and (3.31) as . The edge condition (1.5) for the Dirichlet case implies that and . Using the well-known fact that for any function behaving like as , its Laplace transform , as defined by (3.16), behaves like as , and noting that as , we can show that as , , and within , while within . We can also determine that and as . This means that all parts of equations (3.30) and (3.31) are decaying as in the appropriate half plane. Construct the two functions,
[TABLE]
Both and are entire and decaying at infinity, therefore Liouville’s theorem can be applied to show that . It is hence possible to determine and ,
[TABLE]
Equations (3.33) and (3.34) can be substituted into (3.24) to obtain,
[TABLE]
which is the exact spectral function (2.13) obtained using the S-M technique.
3.4 Neumann boundary condition
The Neumann problem is solved in a similar way to the Dirichlet problem. The transformed Neumann BCs are , leading to a simplification of (3.19). Again, in the resulting form, the W-H technique cannot be applied directly and we need the useful mapping of Section 3.1 together with the factorisation of and given in the same section. This leads to
[TABLE]
The left (resp. right) sides of (3.36) are meromorphic in (resp. ) but as for the Dirichlet case, there are potential poles on . To counteract this we remove the poles on the right side using the knowledge of the GO component of the solution. Assuming that , the GO components of are,
[TABLE]
As in the Dirichlet case, we only need to remove the pole of at with residue,
[TABLE]
Using this residue, and the fact that and , we can remove this pole from the W-H system (3.36) to get
[TABLE]
In both equations (3.38) and (3.39), the left (resp. right) sides are now analytic in (resp. ). As before, in order to apply Liouville’s theorem, we must determine the behaviour of each part in equations (3.38) and (3.39) as . Using the edge conditions (1.5) for the Neumann case, and the reasoning developed in the Dirichlet case, we can show that, as , and within , while within . We can also determine that and as . As before we can hence construct two decaying entire functions and apply Liouville’s theorem to obtain
[TABLE]
Equations (3.40) and (3.41) can be substituted into (3.24) to get,
[TABLE]
which is the exact spectral function (2.18) obtained using the S-M technique.
Critical analysis
The main disadvantage of this method is that, as our derivation shows, it is not naturally designed to tackle the wedge problem. As a result, we do not directly obtain a usual W-H equation, and have to rely on a sophisticated mapping in order to get back to the usual framework.
The advantage of this section, however, is to show the flexibility of the W-H method, and that it can work, even in non-flat/parallel geometries. As for the S-M
technique it is possible to adapt such method to more complicated cases such as inhomogeneous impedance (Shanin, 1998), skew incidence (Daniele and Lombardi, 2006) or even the penetrable wedge (Daniele and Lombardi, 2011).
Moreover, in usual flat geometries, it is known that the W-H technique can be adapted to handle finite structures. It generally results in matrix W-H problems. This is encouraging in our case since there is a chance of tackling geometries such as the truncated wedge (tip removed) with such method. It is perhaps surprising that this problem can be recast in an analytical continuation problem of the W-H problem type. This may give insight to the solution of a broader class of diffraction problems (Daniele and Zich, 2014).
4 The Kontorovich-Lebedev transform method
The third method to be reviewed relies on the Kontorovich-Lebedev (K-L) transform. Introduced in (Kontorovich and Lebedev, 1939), this transform is useful because the resulting transformed Helmholtz equation is easy to solve. However, the inverse transform is known to have convergence issues, but there are alternative versions involving a convergence factor that can help with this (see e.g. (Jones, 1980)). For any function , define the K-L transform and its inverse (which can have many variations from numerous sources such as (Lebedev, 1965), (Abrahams, 1986) and (Jones, 1986)) as
[TABLE]
where and are the Bessel and Hankel functions of the first kind. The transform is valid if,
[TABLE]
for all and . Alternatively, if the second integral condition is not satisfied because tends to a constant as , then contains a pole at on the integration contour of the inverse transform. This pole is interpreted as,
[TABLE]
If the integrand of the inverse transform fails to converge as , an alternative version with a convergence factor (proposed by Jones (1980)), should be used:
[TABLE]
To adapt this to our problem, we first split the total wave field into its incident and scattered parts , where the scattered part satisfies Helmholtz’s equation (1.1) and two types of BCs,
[TABLE]
For our problem, the K-L transform and the associated inverse are given below,
[TABLE]
where the first integral condition (4.2) is satisfied due to the radiation condition (1.4). The edge condition (1.5) implies that will have a pole at . Using (4.7), we find the transformed boundary data,
[TABLE]
From equation 6.611.5 of (Gradshteyn and Ryzhik, 2014), for , we know that
[TABLE]
and integrating (4.10) with respect to , we obtain
[TABLE]
Now, let and , then use (4.11) (resp. (4.10)) to evaluate (4.8) (resp. (4.9)) explicitly to get
[TABLE]
The advantage of the K-L transform is that if satisfies the following governing equation,
[TABLE]
then satisfies Helmholtz’s equation. For the Dirichlet case, the solution of (4.14) is,
[TABLE]
This means that the exact solution with Dirichlet BCs is
[TABLE]
For the Neumann case, the solution of (4.14) is,
[TABLE]
This means that the exact solution with Neumann BCs is
[TABLE]
While it is difficult to see it by inspection, these integral solutions (4.16) and (4.18) are equivalent to the Sommerfeld integral solutions (2.14) and (2.19). As we discuss in B, the connection between the Sommerfeld inverse formula (2.1) and the K-L inverse transform (4.1) was made in (Malyuzhinets, 1958c) by using the Sommerfeld integral form of Bessel functions. Here, we shall show equivalence by first rewriting the integrals (4.16) and (4.18) as series, then convert the result into Sommerfeld integrals.
The solutions (4.16) and (4.18) are evaluated by deforming the contour to the right and summing the residues of the poles crossed in the process. The double pole at is interpreted using (4.3). We can simplify the result using the Jacobi-Anger expansion of the incident wave to obtain the following series solutions for Dirichlet and Neumann BCs respectively,
[TABLE]
where, as before, . These series solutions can be matched with classical series solutions obtained by Macdonald (1902) (see A for more details).
Finally, we need to transform (4.19) and (4.20) into Sommerfeld integrals. We do this by using the Sommerfeld integral formula for the Bessel function of the first kind,
[TABLE]
and equation 1.461.2 from (Gradshteyn and Ryzhik, 2014)
[TABLE]
which converges if . This means that (4.19) and (4.20) can be written in Sommerfeld integral form as
[TABLE]
where the plus and minus signs denote the Dirichlet and Neumann solutions respectively. Using standard trigonometric identities, it is trivial to show that the square brackets in (4.23) are an alternate form of . Hence the Kontorovich-Lebedev solutions (4.16) and (4.18) match with the Sommerfeld integrals (2.14) and (2.19) respectively.
Critical analysis
The main advantage of such transform is that it is a very natural way to tackle the wedge problem. It hence leads to a constructive proof of the form of the solution in the K-L space, see e.g. (4.15). In addition, it leads easily to a near-field expansion of the solution.
The clear disadvantage of such method lies in the convergence issues of the inverse K-L transform. It does require some regularisation in order to be evaluated numerically, and even in this case, the computation of inverse K-L transform remains cumbersome. To evaluate the far-field it is usually necessary to convert it into a Sommerfeld type integral.
5 Solution analysis and evaluation
We shall now compare the exact integral and series solutions with some GTD approximations and evaluate them for some representative values of . Note that the K-L integrals do not need to be plotted since we have already shown that they are equivalent to the Sommerfeld integrals (2.14) and (2.19). Numerical computation of the Sommerfeld integrals can be slow if because will oscillate rapidly along the Sommerfeld contours.
Another way to evaluate these integrals is to deform the Sommerfeld contours to the steepest descent contours shown on the left side of Figure 6. During this deformation, all poles on the real line segment are crossed. Their contribution, which can be calculated exactly using residues, correspond to the GO component of the field, , leaving behind the diffracted part . The steepest descent contour is repeated twice in opposite directions so is cancelled out and the other two contours are translations of each other. The exponential term does not oscillate along these contours so computation time is significantly reduced, even for large . This means that the S-M integrals (2.14)-(2.19) are equivalent to,
[TABLE]
Since they are translations of each other, we can transform and onto the contour which is illustrated on the right side of Figure 6,
[TABLE]
By the method of steepest descent, satisfies,
[TABLE]
where . This is rewritten in terms of the Gudermannian function ,
[TABLE]
Using the following parametrisation, , and noting these identities,
[TABLE]
the diffracted part is written as a simple integral:
[TABLE]
As stated earlier, this integral will be much faster to evaluate numerically than the S-M integrals (2.14)-(2.19). However, difficulties can arise when is in a small neighbourhood of the GO discontinuities because one of the poles will be very close to the contour of integration, which will cause numerical issues.
5.1 Comparison with simpler problems
In this subsection we will show that the solution is consistent with the simple case when the wedge opens up to form a half-space or closes to make a half-plane. First, we look at the case where to form a half-space problem. The solution is easily obtainable via the method of images,
[TABLE]
Here the upper and lower signs denote the Dirichlet and Neumann solutions respectively. Obviously (5.6) is equal to the GO component so we need to show that the diffracted part (5.2) is identically zero. Expressing in terms of the cotangent, we find that,
[TABLE]
which is identically zero due to the periodicity of cotangent. This implies that , as required.
For another comparison, we look at the case where the wedge closes to form a half-plane. Hence we let and match the S-M integrals (2.14)-(2.19) with the known solution to the half-plane problem in terms of Fresnel integrals,
[TABLE]
where the upper and lower signs correspond to the Dirichlet and Neumann solution respectively, is the reflected wave and is the Fresnel integral defined555Fresnel integrals can be written in many different ways, see for example (Abramowitz and Stegun, 1964; Noble, 1958; Assier and Peake, 2012b) and references therein. by
[TABLE]
Having implies that and,
[TABLE]
Hence, we can rewrite the S-M integrals (2.14)-(2.19) in the following form,
[TABLE]
where
[TABLE]
It is possible to express (5.12) in terms of a Fresnel integral (a procedure to do this can be found in section 5.3 in (Babich et al., 2007)), leading to
[TABLE]
Using this, we recover exactly (5.8) from (5.11), as expected. Now that the solution matches with that of the half-space and half-plane problems, we shall focus on deriving the GTD approximation for non-degenerate wedges.
5.2 Geometrical Theory of Diffraction (GTD)
Keller (1962) defined the Geometrical Theory of Diffraction to be an extension of classic Geometrical Optics including diffraction terms. The GTD approximation is simply an asymptotic approximation of the total wave field as , creating a high-frequency or far-field approximation. To derive the GTD approximation of the case presented here, we continue with the method of steepest descent applied to (5.2) as . Equation (5.2) is of the form , where is a big parameter, and . The latter has a saddle point at and is such that . Since is also not zero, we can apply the method of steepest descent in its simplest form (see e.g. (Bleistein and Handelsman, 2010)) to get
[TABLE]
Hence we can write
[TABLE]
In this GTD approximation, the term,
[TABLE]
is known as the diffraction coefficient. Unfortunately, this GTD approximation is singular for certain values of , for example in the case where , the GTD is invalid at , and , which correspond to the GO discontinuities. This is the main issue with GTD: while it is a much more accurate approximation than the Geometrical Optics, it becomes invalid at the GO discontinuities. The pursuit of an approximation that is uniformly valid for all has led to the improved Uniform Geometrical Theory of Diffraction (Kouyoumjian and Pathak, 1974). We follow section 5.5 in (Babich et al., 2007) to find the uniform GTD approximation (UTD).
Restricting ourselves to the specific case where and , there are only two values where the standard GTD is invalid, and . To produce the uniform approximation, we first construct a function that is a linear combination of and defined by (5.12). The idea is to remove the poles causing the singularities in (5.15) and then use the method of steepest descent. Consider the following,
[TABLE]
where is defined in (5.12). The upper and lower signs denote the Dirichlet and Neumann solutions respectively. The combination of and has effectively removed the poles at and , but the pole at remains for all values of . We use the method of steepest descent to approximate ,
[TABLE]
We rearrange (5.17) and use (5.13) and (5.18) to find the UTD approximation.
[TABLE]
where and are the reflections of the incident wave from the top and bottom face respectively. If we restricted the incident angle to instead, we would need to use the following function,
[TABLE]
where the same method as the first case will produce another UTD approximation for
[TABLE]
These two approximations are uniformly valid for , however the BCs are only satisfied in the limit . Another potential accuracy issue occurs when approaches . This situation corresponds to a transition in the GO field, from a case when only one reflected wave is present to a case when two reflected waves occur. Finally, note that using the asymptotic expansions for large argument for the Fresnel integrals will simplify the above formulas to produce the GTD approximation (5.15) again.
5.3 Graphical comparison of evaluation methods
The exact solution to the perfect wedge problem has been written as a Sommerfeld integral on the usual Sommerfeld contour as in (2.14)-(2.19) or on its steepest descent contour as in (5.5). Both formulations are exact and equivalent, but the latter is much easier to evaluate numerically. We have also presented three different approximations, a truncated infinite series (4.19), a GTD approximation (5.15) and a UTD approximation (5.19) or (5.21). In this subsection, we will plot the exact solution and each of the approximations and compare their accuracy and computational speed. For the series solutions we shall truncate at 100 terms, which is enough for the wavenumbers considered here.
In Figure 7, we will consider the wedge defined by for zero incidence angle, . This corresponds to a case where the GO part of the field exhibits two reflected waves. In Figure 8, we consider the same wedge, but with an incident angle , corresponding to a GO field with a single reflected wave. In both cases, we will plot the real part of the total field against for different values of and different BCs. In both figures, the thick plain line represents the exact Sommerfeld solution (SI/SDC), the thick dashed line is the truncated series approximation, the dotted line and the thin line represent the UTD and GTD approximations respectively.
In both Figures 7 and 8, we confirm that,
- •
The series solution is very accurate despite the truncation. However if we want to consider larger values of , more terms will be required to remain accurate, which will slow down its computation.
- •
The GTD approximation has the least overall accuracy and becomes invalid when is close to any GO discontinuities , and . It does however satisfy the correct BCs.
- •
The UTD approximation is a clear improvement to the standard GTD approximation away from the boundaries, in particular it does not have any singularities, but fails to satisfy the BCs.
Both the GTD and UTD approximations appear to improve their accuracy as gets larger. To show this, we take the Dirichlet case with and look at the quantities and against for . Figure 9 (left) illustrates the GTD error and shows that it is a good approximation, provided that is large enough and is not too close to one of the singular angles and (which are indicated by a thin vertical dashed line). In Figure 9 (right), it is clear that the UTD error at the boundary decreases significantly as increases, rendering it a very good approximation everywhere if is large enough. We also reconfirm that the UTD approximation is a large improvement in comparison to the GTD approximation.
Finally, for completeness, we replicate some plots from existing literature using the UTD approximation. Specifically, we replicate the first and last plots of figure 5 in (Hacivelioglu et al., 2011) which is a comparison of an alternate definition for (5.5), the series solution with 100 terms and a similar UTD approximation. In order to replicate these plots, we need to adapt to their geometric configuration by making the substitutions and . We use (5.21) with and . Figure 10 (left) is the Dirichlet case with . Figure 10 (right) is the Neumann case with .
6 Alternative methods
Sections 2, 3 and 4 cover methods that are most commonly used in diffraction theory. In this section, we will briefly present three alternative methods that have been tailored to tackle the perfect wedge problem.
6.1 Embedding Formula technique
The first method to be reviewed is based on the idea of embedding. This idea is relatively new in diffraction theory (Williams, 1982), and has mainly been used for planar cracks and slits, as well as parallel combinations of these (Gautesen, 1983; Martin and Wickham, 1983; Biggs, 2001, 2002). Though, recently, in (Craster and Shanin, 2005) it was adapted to wedges with a rational angle. We will here attempt to summarise this method and consider again our wedge region characterised by . We seek the total field satisfying the Helmholtz equation (1.1), subjected to Dirichlet (1.2) or Neumann (1.3) BCs, as well as radiation and edge conditions (1.4) and (1.5) for a plane wave incidence , with incident angle . The aim of the method is to recover the diffraction coefficient of the diffracted field .
The diffraction coefficient
Using classical separation of variables in the polar coordinates and the edge conditions, it can be shown that the total field has an eigenfunction expansion of the form
[TABLE]
where and is a product of Bessel functions and some trigonometric functions of satisfying the BCs666The multiplicative factor is just here to compensate the near-field behaviour of the Bessel functions, and, doing so, somehow normalise the expansion.. In the Dirichlet case, the term in the sum is equal to zero. Note that using the series results (4.19) and (4.20) of Macdonald type, we can recover exactly, but we will not use this here. The aim is to determine the diffraction coefficient , already defined in (5.16), that is such that
[TABLE]
The edge Green’s functions
In order to do this, as is customary with embedding, we need to introduce an auxiliary problem777Here the auxiliary problems will be constructed from point sources. However, another type of embedding formulae can be obtained with plane wave auxiliary problems, see (Biggs, 2006) for example.. In fact here, we will introduce infinitely many of them. Let , and consider the function that is the tailored Green’s function (i.e. that satisfies the BCs) for the Helmholtz equation resulting from point sources given by , where for Dirichlet BC and for Neumann BC. The strength of each source is given by , as illustrated in Figure 11.
The th edge Green’s function is then defined by
[TABLE]
The near field behaviour of the edge Green’s function can be studied by considering for fixed , close to the wedge edge. In that vicinity, we can scale the space variables to show that behaves locally like , which is the exact same Green’s function, but for Laplace’s equation instead of Helmholtz. Using the method of images in a half-space, and the mapping , it is possible to find explicitly888Note that in (Craster and Shanin, 2005), only the Dirichlet formula is given, and is slightly different from this one (the factors are missing), which we think is a typographical error. as
[TABLE]
where , being the angle measured from the bottom face of the wedge, so that we have . Looking at the leading order of as , using the fact that , we get
[TABLE]
which by construction, is also the local behaviour of near the wedge edge. Note that the edge Green’s function is singular on the wedge edge and does not satisfy the edge condition, we say that it is oversingular. It does however satisfy the Helmholtz equation everywhere outside the wedge. This leads to the exact representation of as:
[TABLE]
since it is clear that the above expression has the right type of singularity, and satisfies the boundary and radiation conditions, as well as the Helmholtz equation.
It is also natural to define the directivity for each edge Green’s function by
[TABLE]
and using the asymptotic behaviour of the Hankel function for large argument, (6.9) and (6.10) imply that
[TABLE]
It is important to note the main difference between the directivities and the diffraction coefficient : the former only depends on one angular variable, while the latter depends on two. Remarkably, using the reciprocity principle, it is possible to relate the far-field of the edge Green’s functions to the near-field of each components of the eigenfunction expansion (6.1) as follows:
[TABLE]
The operator
As mentioned above, this method can only be applied to rational angles999Shanin and Craster (2010) have extended this work by considering a pseudo-differential operator instead of the differential operator . Note that for an integer , reduces to for some constant , which establishes the link with the theory developed here. This new operator can however be used when to produce an embedding formula valid for wedges with non-rational angles, though it cannot be used for polygons., so let us set for some positive integers and . Now define, the operator as follows:
[TABLE]
where is the th Chebyshev polynomial, and it is understood that for some integer , . From now on, for brevity, we will focus solely on the Dirichlet case. It is relatively easy to show that for every , satisfies the Helmholtz equation, the correct boundary conditions and the radiation condition, and that . It is also possible to prove (though it is more difficult) that
[TABLE]
We refer to (Craster and Shanin, 2005) for the details of the proof, but it relies on a careful analysis of the near-field and far-field behaviour of . It also uses the identity , which explains how enters the scene.
Embedding formula
Note now that the behaviour of each term in (6.15) reminds of that of the th edge Green’s function (see (6.6)). This motivates the introduction of the auxiliary function
[TABLE]
By construction, satisfies the edge condition, and it is also clear from what has been done above, that it satisfies the Helmholtz equation, the boundary and the radiation conditions. Hence, by uniqueness, we conclude that , and we obtain the weak form of the embedding formula
[TABLE]
valid everywhere, that relates the total field to the edge Green’s functions. Focusing now on the far-field, (6.16) makes it possible to express the diffraction coefficient in terms of the directivities of some of the edge Green’s functions, as summarised in the equation below:
[TABLE]
The formula (6.17) is the main result of (Craster and Shanin, 2005) and is referred to as the Embedding formula. It is remarkable in the sense that it allows to express the diffraction coefficient, depending on two angular variables, in terms of a sum of products of simpler directivities depending on one angular variable only. Moreover, in that case, thanks to (6.13), we know the directivities exactly and we can then recover a new analytical expression for the diffraction coefficient. For a given rational angle, it is possible to show that it is indeed equal to that given in (5.16).
Critical analysis
The concept of embedding is very general in diffraction theory, which makes this method very adaptable to all kinds of geometries such as slits, wedges, plane sectors, cubes (Skelton et al., 2010) and curved geometries (Moran et al., 2016). In that respect, instead of being seen as a method, one can consider the embedding structure as an inherent property of diffraction problems.
Its main advantage is that once derived explicitly, the embedding formula of a given diffraction problem allows one to obtain a very efficient way of computing the diffraction coefficient resulting from a incident plane wave for all observation and incident angles.
Even though the weak embedding formulae of the type (6.16) are valid everywhere, the power of embedding formulae only becomes apparent in the far-field, where it can be written in its strong form (see (6.17) for the present case). In that sense, such formula is not particularly helpful to shed some light on the near field behaviour of diffraction problems. Though because of this emphasis on the far-field, one can consider structures with multiple diffracting parts such as polygons for example.
6.2 Random Walk method
This method, developed in a series of papers (Budaev and Bogy, 2001, 2002a, 2002b), and applied to the wedge problem in (Budaev and Bogy, 2003), is based on the Feynman-Kac formula (see e.g. (Feynman, 1948), (Kac, 1949) and (Freidlin, 1985)). This formula implies, in particular, that the solution of a deterministic PDE on a domain with Dirichlet condition on the boundary
[TABLE]
with real-valued coefficients , and , can be written as
[TABLE]
where E represents the mean operator, and are random motions governed by the two coupled stochastic differential equations (SDE) with drift coefficient and diffusion coefficient
[TABLE]
with initial conditions (ICs) and , where are Brownian motions (also known as Wiener processes)101010We do not intend to insist on the rigorous mathematical definitions of these objects here, however we refer the interested reader to general textbooks on the topic, such as (Voss, 2013, Chapter 6) for example, where Brownian motions, SDEs (and their resolution via the Euler-Maruyama scheme) and Itô calculus are introduced.. The exit time is the time when each computation should be stopped and it corresponds to the first time such that .
If the coefficients in (6.18) and (6.20) are complex-valued (which as we will see will be the case for the problem at hand), then the Feynman-Kac representation is still valid, but it becomes difficult to determine and define the exit time . In fact if the coefficients are complex, then so will be the random motions , and since the points of belong to , we cannot easily characterise the fact that hits this boundary. This can be addressed by considering the “continuation” of the boundary into a manifold of real dimension 2 within the space and by multiplying (6.18) by , where is a complex-valued function. For a suitable function , it becomes possible to define an exit time , and the solution to the PDE is given by
[TABLE]
where are random motions governed by the two coupled stochastic differential equations (SDE)
[TABLE]
with ICs and .
In order to fit within this framework, for the wedge problem in (Budaev and Bogy, 2003), the authors aim to solve the Helmholtz equation (1.1), subject to the radiation condition and to Dirichlet BCs of the type . They seek a solution of the form for some unknown functions and . The Helmholtz equation becomes
[TABLE]
If the solution is in the Liouville form, we choose . In this case, automatically satisfies the eikonal equation and, after multiplication by , (6.22) becomes
[TABLE]
and we can write the BCs . This fits exactly within the realm of (6.18), but with complex coefficients. It is shown in (Budaev and Bogy, 2003) that a suitable choice of the function is . The manifold extending the boundary is chosen as , where
[TABLE]
In this particular case, the two SDEs to consider become
[TABLE]
with ICs and and exit time defined such that , which now makes sense since the random process is now real for all times. Using (6.21), we can hence write the solution111111Note that in (Budaev and Bogy, 2003), in their equivalent of the second part of (6.24) (their equation (26)), the argument of the exponential is . We believe it to be a typographical error. as
[TABLE]
where the second part of (6.24) is derived from the first using Itô calculus. Note that for this to be valid, should be chosen such that it can be analytically continued for . The second arguments in (6.24) do not pose any problem since by definition . The two SDEs (6.23) are reasonably straightforward to solve numerically (see Figure 12, left) using Euler-Maruyama with time step . If the solution we are trying to find is continuous everywhere, the solution (6.24) can be implemented and works well. To illustrate this point we use the same example as in (Budaev and Bogy, 2003) and apply this method to reproduce the function , which is well known to satisfy the Helmholtz equation and the radiation condition. In order to do so we tailored the BCs to be and plotted an illustration of the result in Figure 12.
If the solution we are seeking has some discontinuities, then the method should be adapted slightly. We are interested in this since what we want to compute is the diffracted field resulting from an incident plane wave with incident angle , which satisfies homogeneous Dirichlet BCs and the radiation condition. In what follows, we choose such that both wedge faces are illuminated. As shown in Section 5, the field has GO discontinuities121212Note that in (Budaev and Bogy, 2003), the convention to choose the index of or is different, but we have made that choice in order to be consistent with the rest of the review. at and , and the knowledge of the field allows us to derive the following jump conditions across :
[TABLE]
where is defined such that , and the bracket . Using these jump conditions, it can be shown that (6.24) can be rewritten as
[TABLE]
where for each simulation, the represent the times of crossings between and the discontinuous lines . If (resp. ) is crossed, then (resp. 1). If the crossing is from above (resp. below), then (resp. ). As illustrated in Figure 13 (left), many such crossings can occur before the exit time is reached. The method has been implemented for a wedge characterised by , and the results, obtained for 2000 realisations (simulated by Euler-Maruyama with time step ), are shown at an observation point , for . Note that if the method was described in (Budaev and Bogy, 2003), it was only implemented for a half-plane, and not for a wedge. Though, as predicted in (Budaev and Bogy, 2003) the error is of the order of 0.01 and the method works well131313 Note that in (Budaev and Bogy, 2003), there is a factor in front of E in the formulae (6.25). This was a typographical error. .
Critical analysis
This method has the advantage of being very adaptable to all sorts of geometries since it is based on the Feynman-Kac theorem that is a very general result (both in terms of geometry and in terms of equation). This adaptability is confirmed by the fact that it has been used in the context of cones (Budaev and Bogy, 2003), quarter-plane (Budaev and Bogy, 2005) and other geometries.
It has to be said however, that for the Helmholtz equation, the PDE and SDE coefficients become complex. This renders the determination of the end time rather more complicated than the real coefficient case. It necessitates to find a convenient complex coefficient to multiply our equations by, and also to find a way of somehow extending the real geometries in a higher dimension complex space.
This method can also become very computational very quickly. Indeed, if one would like for example to recreate a heat map similar to those presented in Figure 15, one would need about 2000 simulations of the SDE system per point, which for a good resolution may lead to a very long computational time.
Another comment that can be made about this method, is that it stands out from all the other methods presented here in terms of the type of mathematics used. This can be considered as an advantage for researchers open to exploring many areas of mathematics, though, this also means that for the usual specialists in diffraction theory, this may result in a steep learning curve.
6.3 The method of functionally-invariant solutions
The third and final alternative method to be reviewed is also known as the Sobolev-Smirnov method. Some recent publications using this method include (Komech et al., 2015; Babich, 2015). The former studies wedge diffraction with a number of different combinations of Dirichlet and Neumann BCs, while the latter studies the impedance wedge problem.
The idea behind this method is to identify the time-harmonic problem with an elementary time-dependent problem where the incident plane wave is replaced with a Heaviside step function such that no diffraction occurs before the time . This means that the solution to this elementary problem (call it ) satisfies the following conditions,
- •
The governing equation is the linear wave equation .
- •
Dirichlet or Neumann BCs at .
- •
can be linearly decomposed into incident and scattered parts, , where .
For simplicity, we shall restrict values of the incident angle and the wedge angle such that, . This restriction means that for , the incident wave does not reach the wedge until when it first touches the wedge at its corner. For , the incident wave has passed the wedge corner and reflected and diffracted waves have appeared. Figure 14 describes this configuration and gives known values of outside the diffraction circle which are found by Geometrical Optics.
The radius of the diffraction circle is and we call the unknown solution inside, i.e. within the diffraction disc, . We need to look at a particular class of solutions to the wave equation and express in terms of a complex variable. Noting that a real solution is required, we write
[TABLE]
ensuring that the wave equation is automatically satisfied. Sections 52-53 in (Smirnov, 1964) gives a detailed explanation as to why this is the case. Note also that within the diffraction disc (i.e. ), the pre-exponential factor in (6.26) is real and positive, varying from [math] to . As a result, maps the diffraction disc onto the unit disc . Therefore, we need to find a function that is analytic inside the unit disc and has the following boundary values,
[TABLE]
where the arc and line numbers are given in Figure 14.
We will now use a conformal mapping to transform this problem into a Riemann-Hilbert problem. In order to do so, we define , where as before , which transforms the indented unit disc of the plane onto the unit upper-half semi-disc of the plane, as illustrated in Figure 14. The branch of the root is defined such that the cut is on the negative real axis and . Let , then we analytically continue into the unit lower-half semi-disc by Schwarz reflection principle (see Figure 14) using anti-symmetry (Dirichlet case) or symmetry (Neumann case). Let and then has the following boundary values for the Dirichlet and Neumann cases,
[TABLE]
The method to solve these two Riemann-Hilbert problems is detailed in section 54 in (Smirnov, 1964). The respective solutions to (6.29) and (6.32) are,
[TABLE]
where the used logarithm has a branch point at with a branch cut along the positive real axis. Using this solution, it is easy to recover the physical solution inside the diffraction disc, and hence the whole solution . We will now see that using a simple Fourier transform, we can recover the sought-after time-harmonic problem from this solution . Consider the evaluation to the following integral, assuming that has a small positive imaginary increment so that is zero,
[TABLE]
With this in mind, we can determine the total field from by using a similar integral,
[TABLE]
and thus, we have found the solution to the time-harmonic problem.
Critical analysis
All of the methods presented so far were tailored to the time-harmonic problem, this means that if one is interested in a time-dependent problem, using these methods would involve taking the inverse Fourier transform of our time-harmonic solutions, which can be expensive computationally. This present method however is tailored to the time dependent problem, which is great if one is interested in the tracking of wave fronts in time for example. It means however that if one is interested in the time-harmonic problem, one would have to take the Fourier transform in time of the solution, as per (6.36), which can prove quite expensive numerically. This method, though in essence designed for the wedge geometry, has been shown to be adaptable to various BCs. One can refer to (Babich, 2015) for example for the case of impedance BCs.
7 Final plots and conclusions
In this review article, we have discussed six different methods that have been applied to the problem of diffraction by wedges with perfect Dirichlet or Neumann boundary conditions. The three main methods discussed were the Sommerfeld-Malyuzhinets technique, the Wiener-Hopf technique and the Kontorovich-Lebedev transform technique. The three alternative methods reviewed were the embedding formula, the random walk method and the method of functionally-invariant solutions (Sobolev-Smirnov). We also looked at two approximation methods, the Geometrical Theory of Diffraction and the Uniform Geometrical Theory of Diffraction.
This list is by no means exhaustive and we should also mention Budaev’s method for elastic wedge scattering (Budaev and Bogy, 1998), the Physical Theory of Diffraction (Ufimtsev, 2014) and an interesting method called the Wiener–Hopf–Hankel formulation (Teixeira, 1991; Castro and Kapanadze, 2010). We note that (Israilov, 2013) could also be applied to the wedge geometry.
We evaluated numerically the exact solution and the associated approximations (series, GTD, UTD) for several configurations and studied their relative performances. We found that the best way to evaluate the exact solution was to consider the integral defined on the steepest descent contour. As regard to the approximations, the truncated series solutions performs very well with low wavenumbers, and we found that while the UTD approximation takes longer to compute, it is a better approximation compared with the GTD because it is uniformly valid and more accurate at lower values of . It has however two main disadvantages, the inaccuracy at the wedge boundary, and also the fact that (5.19) and (5.21) are not continuous across the Geometrical Optics limit .
As emphasised in the critical analysis of each section, the use of the six techniques in this review is not limited to the perfect wedge problem. Examples of extensions include for example impedance wedges (Malyuzhinets, 1958a; Babich, 2015), penetrable wedges (Rawlins, 1999; Lyalinov, 1999; Daniele and Lombardi, 2011) and quarter-plane diffraction (Shanin, 2005; Assier and Peake, 2012a; Budaev and Bogy, 2005; Lyalinov, 2013).
To conclude this review, using the UTD approximation, we produce some density plots of the real part of the diffracted field and the total field in Figure 15 for a wedge defined by and two incident waves, and and a wavenumber . As expected, we see clear discontinuities in the diffracted wave, which are due to GO discontinuities. For the total field, as expected, we see the GO behaviour in the relevant regions, the boundary conditions, and a decaying diffracted field.
Acknowledgements
Nethercote acknowledges the financial support from EPSRC (DTA studentship). Assier acknowledges that part of this work was supported by EPSRC (EP/N013719/1). Abrahams acknowledges the support by UKRI under grant no EP/K032208/1. The authors thank the anonymous reviewers for their help in improving this work. Assier would also like to thank A.V. Shanin for fruitful discussions, especially on the topic of B, and B. Budaev for being so responsive when queried about detailed aspects of the random walk method.
Appendix A Macdonald’s series solution
In this section we shall briefly discuss the separation of variables method applied by Macdonald (1902) to the wedge problem with line source incidence. After this, departing slightly from Macdonald’s approach, we shall use a limiting procedure in order to recover the series solutions (4.19) and (4.20) to the plane wave incidence problem.
The wedge problem forced by a line source of strength with polar coordinates , has the following governing equation,
[TABLE]
where is the Dirac delta function. The total field, , is decomposed into incident and scattered parts where the incident wave is given by
[TABLE]
and is subjected to BCs, (1.2) or (1.3). Considering the ansatz , using separation of variables and applying the BCs, we obtain the following series solutions:
[TABLE]
Because of the source location at , and the need to satisfy both the edge and the radiation conditions (satisfied by the Bessel and Hankel functions respectively), we pose
[TABLE]
To ensure continuity across , we require and . We can determine the coefficients and by deriving and applying a jump condition across .
In the Dirichlet case, substitute (A.3) into (A.1), and multiply the resulting equation by . Integrating w.r.t. from to , and using the orthogonality of sine, we obtain
[TABLE]
Now integrating (A.6) from to and taking the limit leads to the jump condition
[TABLE]
Lastly, we use (A.5) and the Wronskian result to determine that . Hence, the series solution with line source incidence and Dirichlet BCs is
[TABLE]
where and . This agrees with Macdonald’s solution141414Note that Macdonald (1902) uses the alternate time factor ..
For the Neumann case, the coefficients are found by the same method using the orthogonality relation for cosine, leading to
[TABLE]
where and for .
To recover the plane wave solution, we send the source and its strength to infinity in a way that ensures that (as defined in (A.2)) behaves like as . This can be done by choosing and leads to
[TABLE]
Hence, for plane wave forcing with Dirichlet or Neumann BCs respectively, the series solutions are
[TABLE]
which matches perfectly with (4.19) and (4.20) as required. Note that these exact series solutions have a natural embedding structure (see Section 6.1) in the sense that they are simply sums of products of functions of one variable only.
Appendix B A link between the spectral function and Green’s integral operator
B.1 Preliminary definitions and Green’s function representation
Let us introduce the generic plane wave function by
[TABLE]
It is important to note that for any , is exponentially decaying as as long as , where both and are understood as open sets (do not contain their boundaries) and are illustrated in Figure 16.
Using the notations of Section 5, the total, diffracted and geometrical optics fields are denoted , and . For the exterior wedge, consists of an incident wave and one or two reflected waves and can hence be written in the form where is either zero or a given constant. Typically in our problem the incident wave corresponds to , and the reflected waves to either or or both depending on how many reflections we have.
For a given function , let us now introduce the Green integral operator defined by
[TABLE]
In the contest of this review, . Moreover, using standard integration, one can show that , and hence we can write
[TABLE]
which implies in particular that each of the are simple poles of with residue . The objective of this appendix is to find a connection between the spectral function and . In order to do that, we shall make use of the theory of Green’s functions as follows.
Let us pick a point , and pick two angles and (the subscripts and stand for above and below) chosen such that , and a radius . Now consider the domain to be the corresponding sector described in Figure 17. Let us further assume that is oriented anti-clockwise, and that the normals to are chosen to be outgoing.
Let be a short notation for , the free-space Green’s function for the Helmholtz equation resulting from a point source at . Using the respective governing equations of and , and the divergence theorem, we have
[TABLE]
Hence, using that on , , on , and on , , we get
[TABLE]
Using the Hankel representation of , its far-field asymptotics, and the method of steepest descent, we can show that the only part of the far-field leading to any contribution of the arc integral as is an incident plane wave coming from within the sector. More precisely, if ,
[TABLE]
All other components (reflected waves, diffracted field, incident waves from outside the sector) can be shown to have zero contribution. Hence, taking the limit as in (B.4), we get
[TABLE]
where if and zero otherwise. Hence the knowledge of and on oblique lines of constant is important. At this stage, it is important to realise, at least informally, that if we could write them in terms of somehow, then we have a chance to link and the Green integral operator.
B.2 Green’s functions on oblique lines
Before finding formulae for , we will focus on the Green’s function corresponding to a point source at the origin. First of all, it is well known (see e.g. Kythe (2011)) that . Moreover, the Hankel function has the following integral representation151515See, e.g. Sommerfeld (2003) eq (6) p19, together with translators’ note 4 on p78, here we use .
[TABLE]
where is described in Figure 6. We will now endeavour to find formulae for valid on an oblique half-space and hence on any line that crosses the axis with an angle say, and lies above (see Figure 18 (left)) or below (see Figure 18 (right)) the origin.
Oblique line above the origin
Let us consider the half-space , the grey area of Figure 18 (left) . Let us start from (B.6) and shift the contour to the contour , where the new contour height is adjusted so that it goes through the origin. Because of the restriction on , we can do that without leaving , where our integrand is analytic and exponentially decaying, and so the value of the integral and its convergence property remain unchanged. We can now perform the substitution to get
[TABLE]
where the contour goes through the point of the real axis, as shown in Figure 19 (left). This formula is valid (and the integral converges exponentially) on any oblique line with angle that lies above the origin.
Note that here is finite, and the integrand is analytic, so we can in principle deform all the contours to , it is important to note that the latter crosses the real axis at and is included (only just!) in . This contour will just be referred to as thereafter, and we get
[TABLE]
Oblique lines below the origin
In a very similar way, we can consider the half-space and shift the contour to a height adjusted passing through the origin. Upon performing the substitution , we obtain an integral over a contour illustrated on Figure 19 (right). Again by analyticity of the integrand, such integral can safely be deformed to the contour that crosses the real axis at and lies within , to get
[TABLE]
Back to
In order to get back to , we just need to replace by and by in (B.8) and (B.9), where and are the polar coordinates centred at . Upon noting that , we find that
[TABLE]
where from now on, the subscript is either or . Since is independent of , we get a similar formula for . In particular, in the configuration of Figure 17, since the oblique line (resp. ) lies above (resp. below) the source and make an angle (resp. ) with the real axis, we have
[TABLE]
B.3 Connection formula between and
Before making use of our results (B.5) and (B.10), we need to make use of some properties of the Green integral operator:
Proposition 2
**
Apart from eventual poles on the real line, as a function of , is analytic for . 2. 2.
If , then . Note that by analytic continuation, this allows to extend the natural domain of analyticity of to .
Now, we can input (B.10) into (B.5), and, since we made sure that , we can exchange the order of integration. Let us furthermore assume that , then the formula can be evaluated at to get
[TABLE]
where an illustration of the contour configuration is displayed in Figure 20 (left). Making use of point of Proposition 2, the integrands of both integrals in (B.11) are actually analytical continuations of each other, and hence we can write
[TABLE]
where is a notation for going in the other direction. Let us consider the contours as angular (we can do that by analytic deformation), as depicted in Figure 20 (right). Let us also introduce a new contour , that is rectangular, with its centre at the origin and oriented anticlockwise, such that its left (resp. right) lateral side coincides with a part of (resp. ), but in the opposite direction.
We can always choose and close enough to zero, such that no poles related to reflected waves exist within . In this case, one can show that the only possible singularity is a pole corresponding to the incident wave, and we have
[TABLE]
This ensures that we can write
[TABLE]
Now, the coinciding lateral parts cancel each other, and the remaining contour is simply (see Figure 2 (left)). Now taking the limit as , or using the fact that is an analytic continuation of by Proposition 2, we get
[TABLE]
Everything that has been done in this subsection can be used to get a similar formula for to get
[TABLE]
Now, comparing (B.14) and (B.15) to equations (2.1) and (3.22), it is clear that we can apply Theorem 1 to find that
[TABLE]
leading to the sought-after formula
[TABLE]
linking the spectral function to the Green’s operator . Note that this formula could have been recovered from what was done at the end of the Wiener-Hopf section (Section 3) in particular it is a consequence of (3.24), but this appendix is showing this link from Green’s identity only. This can also be seen as a constructive way of getting to the form of the Sommerfeld integral. We can also follow the paper (Malyuzhinets, 1958c) to directly link the spectral function with the Kontorovich-Lebedev transform of the scattered wave ,
[TABLE]
Note that the work done in this appendix is very general and can possibly be applied to geometries other than the wedge.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Abrahams (1986) Abrahams, I. D., 1986. Diffraction by a semi-infinite membrane in the presence of a vertical barrier. J. Sound Vib. 111 (2), 191–207.
- 2Abrahams (1987) Abrahams, I. D., 1987. On the Sound Field Generated by Membrane Surface Waves on a Wedge-Shaped Boundary. Proc. R. Soc. A Math. Phys. Eng. Sci. 411, 239–250.
- 3Abramowitz and Stegun (1964) Abramowitz, M., Stegun, I. A., 1964. Handbook of Mathematical Functions. Courier Corporation.
- 4Assier and Peake (2012 a) Assier, R. C., Peake, N., 2012 a. On the diffraction of acoustic waves by a quarter-plane. Wave Motion 49 (1), 64–82.
- 5Assier and Peake (2012 b) Assier, R. C., Peake, N., 2012 b. Precise description of the different far fields encountered in the problem of diffraction of acoustic waves by a quarter-plane. IMA J. Appl. Math. 77 (5), 605–625.
- 6Babich (2015) Babich, V. M., 2015. Solution of the Diffraction Problem of a Plane Wave by an Impedance Wedge (Non-Stationary Case, Smirnov-Sobolev Method). Russ. J. Math. Phys. 22 (2), 145–152.
- 7Babich et al. (2007) Babich, V. M., Lyalinov, M. A., Grikurov, V. E., 2007. Diffraction Theory: The Sommerfeld-Malyuzhinets Technique. Alpha Science.
- 8Bender and Orszag (1999) Bender, C. M., Orszag, S. A., 1999. Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory. Springer-Verlag, New York.
