Reflections from the open end of a pipe

I made a video about why sound reflects from the open end of a pipe. As usual, this turned out to be more complex than I expected, and the video follows my journey through the physics. This post contains some extra technical notes about details that I glossed over in the video.

Derivation of wave equation for sound

This is based on the derivation in Volume I Chapter 47 of the Feynman lectures; more detail can be found there. I changed some of the symbols to make it look less daunting, for example Feynman uses χ for displacement and I used s.

Simulation of wave equation in 2D

If you look closely at the simulation at 03:14, you will notice that there is a ‘wake’ behind the propagating wavefront: even though the disturbing stimulus is only a short pulse, the pressure decays slowly after the wavefront passes, and the inside of the circle does not return to perfect black. This is a property of wave propagation in 2D which does not happen in 3D (in 3D the wake cancels, and the interior pressure would be zero once the wave has passed). In fact, because of this, the Huygens-Fresnel principle does not strictly hold in 2D, which might come as a shock since you have probably always seen it depicted in 2D. Fortunately the wake phenomenon is subtle enough that 2D problems are still pedagogically useful.

Separation of variables

At 05:35 I start from the assumption of a monochromatic (single-frequency) source at frequency ω that has been going for all time, and thus separate the pressure into a static spatial solution and an oscillation in time:

$\begin{equation*} p(x,y,t) = p(x,y) e^{-i\omega t} \end{equation*}$

I did not state it in the video, but it is useful to know that if p(x,y,t) satisfies the wave equation, then now p(x,y) satisfies the Helmholtz equation:

$\begin{equation*} \nabla^2 p = -k^2 p \end{equation}$

where k is the spatial frequency (=ω/c).

The monochromatic source is just a temporary mathematical tool and not at all limiting. We can describe any other source function as a sum of monochromatic frequency components. And the wave equation is linear, so if we know the solution for each individual frequency, the overall solution is then simply a sum of these solutions.

Decomposition of plane waves

At 06:07, I state that “regardless of what the spatial solution looks like, we can always decompose it into a sum of plane waves using […] a Fourier transform”.

Strictly speaking, this is mathematically true for any function regardless of whether it is a physically meaningful solution of this problem. That is, we can decompose any (well-behaved) 2D function into plane waves by taking a 2D Fourier transform. The problem is, this decomposes the function into plane waves of various spatial frequencies. Most of the waves in this basis do not individually satisfy the above Helmholtz equation (which is only satisfied by waves of spatial frequency k), and that limits the usefulness of this decomposition.

In this case we know that p(x,y) satisfies the Helmholtz equation, and so only waves of overall spatial frequency k are allowed. How do we decompose the solution in this basis? Usually we use what is called the angular spectrum method, which goes as follows for 2D. We choose a line that we know the function on (any line) and take the 1D Fourier transform on that line. Let’s say we take the line x=0 and come up with a spectrum of spatial frequencies k_y. Then for each y-direction spatial frequency k_y there is one corresponding plane wave that has overall spatial frequency k, it has¹:

$\begin{equation*} k_x = \sqrt{k^2 - k_y^2} \end{equation*}$

What happens when k_y > k? We have an imaginary k_x. Since plane waves are of the form:

$\begin{equation*} p(x,y) = Ae^{ik_x x + ik_y y} \end{equation*}$

an imaginary k_x is equivalent to:

$\begin{equation*} p(x,y) = Ae^{-|k_x|x + i|k_y| y} \end{equation*}$

There is a decaying exponential in x instead of a sinusoidal oscillation. It turns out that such waves are still solutions of the Helmholtz equation. They are called evanescent waves. The decaying components perpendicular to the boundary represent ‘unnatural’ forcing conditions at the boundary that cannot be represented as plane waves, and indeed these are precisely the waves that do not radiate.

In case you are wondering: this set of waves (plane waves and evanescent waves of overall frequency k) forms a complete basis to represent solutions of the Helmholtz equation in free space with boundaries. But it is not a complete basis for any function on the 2D space, which would require the full 2D Fourier transform.

The pressure distribution on x=0

At 06:34, I use a rectangular function to approximate the pressure distribution on the line that crosses the aperture (opening). There are a few important comments I need to make about this.

Firstly, a strictly rectangular function is not physically possible. I already said earlier in the video that pressure cannot have discontinuities except at boundaries. Nonetheless, I use this as an example of a spatially constrained function; it is more intuitive to talk about the width of a rectangular function than one with more gentle roll-off.

I also assume that pressure is approximately zero at the wall. While this is presumably true far from the aperture, there is nothing that would force this to hold close to the aperture.

A better (but more complex) approach, then, is to fix the horizontal fluid velocity at the x=0 line instead of the pressure. We do know that the horizontal velocity is zero at the wall, since fluid particles cannot go through the wall. And within the aperture, the horizontal velocity is close to constant for low frequencies².

For each plane wave component that we obtain from a Fourier transform, we have the pressure-velocity relationship³:

$\begin{equation*} \tilde{p} = \rho_0 c \tilde{v} \end{equation*}$

v here is in the direction of propagation of the plane wave. Considering horizontal velocity only, we have v_x = v cos θ and thus:

$\begin{equation*} \tilde{p} = \frac{\rho_0 c}{cos \theta} \tilde{v_x} \end{equation*}$

So if we know the v_x distribution on the line x=0, we take its Fourier transform to get $\tilde{v_x}$ and from the above equation we immediately have $\tilde{p}$ , the pressure Fourier transform. Constants aside, we conclude that the simple sinc function shown in the video should be replaced by one that is corrected by a factor 1/cosθ.

But this causes a bit of a catastrophe as θ approaches 90 degrees. If $\tilde{v_x}$ is non-zero in this
region of the spectrum (which it may be), the plane wave amplitude diverges to infinity.

Actually, this is perfectly reasonable physically. Near-vertical plane waves are going in the near-vertical direction, and contribute negligibly to horizontal velocity. As we approach 90°, the amplitude needs to be ever-larger to drive a given horizontal velocity.

Fortunately, what we care about physically is not the individual plane waves (which are a mathematical decomposition tool), but the power radiated over a given angle. When we work out the power, we find that there is a cos θ in the numerator (see below). This cancels the cos θ in the denominator.

Power

Here I will derive the actual power radiated over a given range of angles in 2D.

Assume we use the following Fourier transform convention:

$\begin{align*} \tilde{p}(k_y) &= \int_{-\infty}^{\infty} p(0,y) e^{-ik_y y} \,dy \\ p(0,y) &= \frac{1}{2\pi} \int_{-\infty}^{\infty} \tilde{p}(k_y) e^{ik_y y} \,dk_y \end{align*}$

We can extend this 1D Fourier spectrum to plane waves over full (x,y) by writing:

$\begin{align*} p(x,y) &= \frac{1}{2\pi} \int_{-\infty}^{\infty} \tilde{p}(k_y) e^{i(k_x x + k_y y)} \,dk_y \\ & \quad\quad\quad\text{where } k_x = \sqrt{k^2-k_y^2} \text{ as before} \end{align*}$

Move to polar co-ordinates (r,θ) where x = r cos θ and y = r sin θ:

$\begin{align*} p(r,\theta) &= \frac{1}{2\pi} \int_{-\infty}^{\infty} \tilde{p}(k_y) e^{i(k_x r \cos \theta + k_y r \sin \theta)} \,dk_y \\ &= \frac{1}{2\pi} \int_{-\infty}^{\infty} \tilde{p}(k_y) e^{i\Phi} \,dk_y \\ & \quad\quad\quad\text{where } \Phi = r(k_x \cos \theta + k_y \sin \theta) \end{align*}$

(Φ represents the phase at (r,θ) and is also a function of k_y.)

Using the stationary phase approximation, we can reduce this to:

$\begin{align*} p(r,\theta) &\approx \frac{1}{2\pi} \, \tilde{p}(k_y^{*}) \sqrt{\frac{2\pi}{|\Phi''(k_y^{*})|}} e^{i(\Phi(k_y^{*}) + ...)} \\ & \quad\quad\quad\text{where } k_y^{*} = k \sin \theta \text{ (stationary point)} \end{align*}$

We are not interested here in the complex exponential (phase), we are only interested in the pressure amplitude:

$\begin{align*} |p(r,\theta)| &\approx \frac{1}{2\pi} |\tilde{p}(k_y^{*})| \sqrt{\frac{2\pi}{|\Phi''(k_y^{*})|}} \\ &= \frac{1}{2\pi} |\tilde{p}(k \sin \theta)| \sqrt{\frac{2\pi k \cos^2 \theta}{r}} \end{align*}$

The time-averaged intensity of the wave (power per unit length in 2D) is related to pressure via:

$\begin{align*} \langle I \rangle &= \frac{|p|^2}{2\rho_0 c} \\ &= \frac{k \cos^2 \theta}{4\pi \rho_0 c r} |\tilde{p}(k \sin \theta)|^2 \end{align*}$

Integrating over small segments of the circle (length r.dθ), the time-averaged power is:

$\begin{align*} P &= \int_{\theta_1}^{\theta_2} \langle I \rangle r \,d\theta \\ &= \frac{k}{4\pi \rho_0 c} \int_{\theta_1}^{\theta_2} |\tilde{p}(k \sin \theta)|^2 \cos^2 \theta \,d\theta \end{align*}$

Again remember this is for 2D, this expression should not be used for 3D. Notice the cos² θ in the integrand; in the previous section we ended up with cos θ in the denominator of $\tilde{p}$ and these will cancel ensuring finite power.

Final comments

At 9:17 I make a throwaway comment about evanescent waves being “not in resonance with the universe”. I think this might be a good way to think about the radiation problem in general. If you think about (say) strings going in every direction that resonate when properly stimulated at a spatial frequency that matches λ=c/f, that might help to demystify why waves can be stimulated in oblique directions, and why high spatial frequencies cannot stimulate waves in any direction. Indeed this is equivalent to calculating the overlap/projection between the local field and the possible plane wave components, which is a classical method of computing radiation, but I do like the idea of drawing strings or something similar as a tool for intuition. I am yet to come up with a good visualization of it however.

Spatial frequencies can be treated like vectors, i.e. a spatial frequency k in some direction can be resolved into a k_x component (= k cos θ) and a k_y component (= k sin θ) with the usual Pythagorean property k_x² + k_y² = k². (Unlike wavelengths which obey the inverse Pythagorean theorem; the wavelengths of a wave in x and y directions are greater than the wavelength in the direction of propagation.) ↩︎
Assume that the left end of the pipe has some sort of piston-like source that drives constant velocity and pressure across the cross-section. For low frequencies where radiation is insignificant, we can assume that this is approximately true right to the end of the pipe. For higher frequencies where radiation is significant, we have coupling of the external field back into the pipe, so the fields across the pipe opening are not constant. The problem can theoretically be solved by decomposing the field at the pipe opening into different orthogonal modes (e.g. one which is constant across the cross-section, one with a single central hump, and so on), calculating the radiation coupling from each mode, and solving simultaneous equations. But this is outside the scope of this post. ↩︎
This can easily be derived from either equation 1+2 or equation 3 in the video, noting that for a plane wave we have δp/δt = –iωp, δp/δx = ikp, etc. p/u = ρ₀c is called the acoustic specific impedance of the medium. Note that acoustic specific impedance gives p/u for individual plane waves, for superpositions of waves p/u is not such a simple function of the medium. ↩︎

Reflections from the open end of a pipe

Derivation of wave equation for sound

Simulation of wave equation in 2D

Separation of variables

Decomposition of plane waves

The pressure distribution on x=0

Power

Final comments

Leave a Reply Cancel reply

Receive new posts by email:

Recent posts:

Categories: