---
jupyter:
  description: 'The Helmholtz resonator with planex2d and NGSolve: the note of a glass bottle with a flat top from an
    axisymmetric acoustic eigenvalue solve, against the classical formula with the end corrections of Norris and Sheng
    and of Levine and Schwinger, for necks from 5 to 100 mm.'
---

# 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.

![The pressure of the resonance in and around a glass resonator, uniform in the bottle and falling along the neck](../_static/gallery/helmholtz-resonator.png){.result}

```python
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
```

```python nbsphinx="hidden"
from docs_figures import MAP, SOFT, Chart, Row


def sampled(p, mesh, box):
    """the mode p for Row, at points of the air inside the view box ((h0, h1), (v0, v1))"""
    def field(face, plot=True, **_):
        (h0, h1), (v0, v1) = box
        gh, gv = np.meshgrid(np.linspace(h0, h1, 141), np.linspace(v0, v1, 281))
        uv = np.c_[gh.ravel(), gv.ravel()]
        uv = uv[face.Distance(uv, signed=True) < -0.05]
        return {"points": [tuple(q) for q in uv], "values": np.asarray(p(mesh(uv[:, 0], uv[:, 1]))).ravel(),
                "label": "pressure, 1 in the bottle"}
    return field
```

## 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):

```python
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
```

```python nbsphinx="hide-input"
air, glass = resonator(50)
Row((air, "A neck of 50", {"fills": [(glass, SOFT)], "view": ((0, 60), (-110, 70))}),
    caption="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:

```python
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")
```

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):

```python nbsphinx="hide-input"
f, mode, mesh = results[50]
box = ((0, 60), (-110, 70))
Row((resonator(50)[0], "The pressure of the resonance", {"field": sampled(mode, mesh, box), "cmap": MAP.reversed(),
                                                          "clim": (0, 1), "view": box, "bar": "right"}))
```

## 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:

```python
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")
```

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):

```python nbsphinx="hide-input"
lengths = np.linspace(5, 100, 96)


def formula(L, d):
    return c / (2 * math.pi) * math.sqrt(A / (V * (L + t / 2 + d * r)))


Chart("neck length in mm", "frequency in Hz",
      dots={"finite elements": (list(necks), [results[L][0] for L in necks])},
      lines={"flanged inside, unflanged outside": (lengths, [formula(L, 0.8216 + 0.6133) for L in lengths]),
             "flanged at both ends": (lengths, [formula(L, 2 * 0.8216) for L in lengths])},
      title="The resonance of the bottle against its neck length")
```

```python nbsphinx="hidden"
import sys
sys.path.insert(0, "..")
from gallery import save_figure
f, mode, mesh = results[50]
air, glass = resonator(50)
axis = ((0, 0), (0, 1))


def both_halves(face, plot=True, **_):
    """the mode for PlotField on the whole section, the half section's values mirrored"""
    gh, gv = np.meshgrid(np.linspace(-135, 135, 541), np.linspace(-110, 70, 361))
    uv = np.c_[gh.ravel(), gv.ravel()]
    uv = uv[face.Distance(uv, signed=True) < -0.05]
    return {"points": [tuple(q) for q in uv], "values": np.asarray(mode(mesh(np.abs(uv[:, 0]), uv[:, 1]))).ravel(),
            "label": "pressure, 1 in the bottle"}


frame = px.QuickSketch((-135, -110)).Rectangle(270, 180).Face()
fig, ax = plt.subplots(figsize=(6, 4))
px.PlotField(both_halves, (air + air.Mirrored(*axis)) * frame, ax=ax, colormap=MAP.reversed(), min=0, max=1)
px.Plot(glass + glass.Mirrored(*axis), ax=ax, face_col="white", line_col="white", arc_col="white")
ax.set_aspect("equal", adjustable="box")
ax.set_xlim(-135, 135)
ax.set_ylim(-110, 70)
save_figure(fig, "helmholtz-resonator", bare=True, bars={1: "pressure, 1 in the bottle"})
```