The Helmholtz resonator

Why does a bottle hum when one blows across its neck? The air in the neck moves as a piston on the spring of the air in the bottle, and the note is that mass on that spring:

\[f = \frac{c}{2\pi}\sqrt{\frac{A}{V L_{eff}}},\]

with \(A\) the neck’s cross-section, \(V\) the bottle’s volume and \(L_{eff}\) the neck’s length plus end corrections for the air that moves with it beyond both ends. Where the neck meets a flat wall, the correction is \(0.8216\,r\) for a neck of radius \(r\), Norris and Sheng’s value for a flanged end, and at a free end it is \(0.6133\,r\), Levine and Schwinger’s value for an unflanged pipe.

Here, NGSolve finds the resonance as the lowest acoustic mode of a glass resonator with a flat top, for necks from 5 to 100 mm long, and the formula is checked against it.

import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import numpy as np
import planex as px
from ngsolve.solvers import PINVIT
from scipy.optimize import brentq

The resonator

The resonator is round, so the problem is axisymmetric: the sketch is the half section with \(h = r \ge 0\), in millimeters. The bottle is 80 wide and 100 high below its flat top, and the neck of radius 10 stands on the top. The wall is a polyline along the glass, and its offset by half of 1 mm on each side is the glass itself, rounded at the lip. The air is a half disc of radius 300 around the resonator, less the glass. Its round edge is named before the Boolean operation, and the axis, the straight edge, too (Figure 1):

R, H, a, t, far = 40.0, 100.0, 10.0, 1.0, 300.0


def resonator(L):
    wall = px.FromPolyline([(0, -H), (R, -H), (R, 0), (a, 0), (a, L)])
    glass = wall.Offset(t / 2)
    glass.edges.maxh = 2.0
    half = px.QuickSketch((0, -far)).Arc(far, 180).Close().Face()
    half.edges.name = "far"
    half.edges(kind="line").name = "axis"
    return half - glass, glass

Figure 1: The air in and around the resonator, the glass in gray

The resonance

The pressure \(p\) of a mode satisfies \(\nabla^2 p + k^2 p = 0\) with \(k = 2\pi f / c\), and the glass reflects the sound, so \(\partial p / \partial n = 0\) on it. Around the axis, the weak form weighs everything by \(r\): \(\int r\,\nabla p \cdot \nabla q = k^2 \int r\,p\,q\). The far edge of the half disc is a sphere of radius 300 with \(p = 0\), which stands for the open air. An empty sphere of radius 300 with \(p = 0\) on it has its lowest mode where \(300\,k = \pi\), at 572 Hz, far above the resonance. PINVIT finds the lowest modes, and resonance returns the lowest frequency and its mode:

c = 343e3


def resonance(L):
    air, glass = resonator(L)
    mesh = ngs.Mesh(px.OCCGeometry(air).GenerateMesh(maxh=50, grading=0.3))
    mesh.Curve(2)
    V = ngs.H1(mesh, order=3, dirichlet="far")
    p, q = V.TnT()
    K = ngs.BilinearForm(ngs.x * ngs.grad(p) * ngs.grad(q) * ngs.dx).Assemble()
    M = ngs.BilinearForm(ngs.x * p * q * ngs.dx).Assemble()
    lam, modes = PINVIT(K.mat, M.mat, pre=K.mat.Inverse(V.FreeDofs(), inverse="umfpack"), num=2, maxit=30,
                        printrates=False)
    mode = ngs.GridFunction(V)
    mode.vec.data = modes[0]
    mode.vec.data = 1 / mode(mesh(1.0, -H / 2)) * mode.vec  # a mode's scale is arbitrary: 1 in the bottle
    return c * math.sqrt(lam[0]) / (2 * math.pi), mode, mesh


necks = (5, 10, 20, 50, 100)
results = {L: resonance(L) for L in necks}
for L, (f, _, _) in results.items():
    print("neck", L, "mm | resonance", round(f, 1), "Hz")

neck 5 mm | resonance 296.6 Hz
neck 10 mm | resonance 266.1 Hz
neck 20 mm | resonance 224.8 Hz
neck 50 mm | resonance 164.0 Hz
neck 100 mm | resonance 122.4 Hz

In the mode of a neck of 50, the pressure is the same everywhere in the bottle, falls along the neck and is nearly zero a little beyond the mouth (Figure 2):

Figure 2: The pressure of the resonance

Against the formula

The air inside the glass is half a millimeter smaller on every side than the wall’s polyline. The neck’s air has the radius \(r = 9.5\) and runs from the bottle’s ceiling to the mouth, \(L + 0.5\) long, and the bottle’s air is 79 wide and 99 high. With the end corrections of a flanged end inside and an unflanged one outside, the formula comes within a few percent of the resonance. Taking the neck’s air as a short tube instead of a rigid piston, \(k \tan(k L_{eff}) = A/V\), the formula comes within a fraction of a percent once the neck is longer than its diameter:

r, V = a - t / 2, math.pi * (R - t / 2) ** 2 * (H - t)
A = math.pi * r**2
for L, (f, _, _) in results.items():
    Leff = L + t / 2 + (0.8216 + 0.6133) * r
    lumped = c / (2 * math.pi) * math.sqrt(A / (V * Leff))
    tube = c / (2 * math.pi) * brentq(lambda k: k * math.tan(k * Leff) - A / V, 1e-6, 1.5 / Leff)
    print("neck", L, "| finite elements", round(f, 1), "| formula", round(lumped, 1),
          "| as a tube", round(tube, 1), "Hz")

neck 5 | finite elements 296.6 | formula 301.7 | as a tube 301.1 Hz
neck 10 | finite elements 266.1 | formula 268.6 | as a tube 268.0 Hz
neck 20 | finite elements 224.8 | formula 225.9 | as a tube 225.1 Hz
neck 50 | finite elements 164.0 | formula 164.8 | as a tube 163.8 Hz
neck 100 | finite elements 122.4 | formula 123.5 | as a tube 122.2 Hz

The shortest neck lies between the formula with a free outer end and the one with two flanged ends. Its mouth is so close to the flat top that the top acts partly as a flange there too, and a flanged end has the longer correction (Figure 3):

Figure 3: The resonance of the bottle against its neck length