Tuning a marimba bar

Why do marimba and vibraphone bars have an arch cut into their underside? A plain bar sounds out of tune with itself, because the first overtone of a free bar rings at 2.756 times its fundamental, not at a musical interval. Thinning the middle of the bar lowers the fundamental much more than the overtone, and an arch of the right depth puts the overtone two octaves up, at four times the fundamental. A vibraphone bar is tuned to these two. Marimba makers also place the next overtone, near ten times the fundamental, with a more shaped cut. The bar hangs on cords through the nodes of its fundamental, where it does not move. On a plain bar, they lie 0.224 of the length from each end.

The example is an aluminum bar like a vibraphone’s, 350 mm long and 13 mm thick, in plane stress. NGSolve and SciPy find its bending modes, which are checked against beam theory for the plain bar, and SciPy’s brentq finds the depth of the arch that tunes the overtone to four times the fundamental.

This page is part of a collection of physics examples using planex2d.

import math

import matplotlib.pyplot as plt
import ngsolve as ngs
import numpy as np
import planex as px
import scipy.sparse as sp
from scipy.optimize import brentq
from scipy.sparse.linalg import eigsh

The bar

The bar is 350 long and 13 thick, in millimeters. It is aluminum, with \(E\) = 70000 N/mm², Poisson’s ratio 0.33 and the density 2.7e-9 t/mm³, so that the frequencies come out in Hz. Its width does not enter, since a bar bending in the plane of its length and thickness has the same frequencies however wide it is. The arch is a disc taken out of the underside. Its circle passes through the bottom edge 0.3 of the length on either side of the middle and reaches depth up in the middle:

L, t = 350, 13
E, nu, rho = 70e3, 0.33, 2.7e-9


def bar(depth, chord=0.6 * L):
    """the bar with an arch of the given depth cut under its middle"""
    face = px.QuickSketch().Rectangle(L, t).Face()
    if depth > 0:
        R = (chord**2 / 4 + depth**2) / (2 * depth)
        face = face - px.Circle((L / 2, depth - R), R).Face()
    return face

Its modes

The bar hangs free, without supports. Its modes are the eigenvectors of its stiffness and its mass, which SciPy’s eigsh finds nearest to a shift below zero. The first three are the bar moving as a rigid body, at zero frequency, and are skipped. The shift is of the size of the eigenvalues, the squares of the circular frequencies. With a shift too close to zero, the rigid-body modes would leave the shifted matrix nearly singular and the frequencies inaccurate:

mu, lam = E / (2 * (1 + nu)), E * nu / (1 - nu**2)


def modes(face, k=8):
    """the bending frequencies in Hz and the modes of the free bar, lowest first"""
    mesh = ngs.Mesh(px.OCCGeometry(face).GenerateMesh(maxh=t / 3))
    V = ngs.VectorH1(mesh, order=3)
    u, v = V.TnT()

    def stress(w):
        strain = ngs.Sym(ngs.Grad(w))
        return 2 * mu * strain + lam * ngs.Trace(strain) * ngs.Id(2)

    a = ngs.BilinearForm(ngs.InnerProduct(stress(u), ngs.Sym(ngs.Grad(v))) * ngs.dx).Assemble()
    m = ngs.BilinearForm(rho * u * v * ngs.dx).Assemble()
    A = sp.csr_matrix(a.mat.CSR(), shape=(V.ndof, V.ndof))
    M = sp.csr_matrix(m.mat.CSR(), shape=(V.ndof, V.ndof))
    w2, vectors = eigsh(A, k=k, M=M, sigma=-1e5, which="LM")
    order = np.argsort(w2)[3:]
    return np.sqrt(w2[order]) / (2 * math.pi), vectors[:, order], mesh, V

The plain bar

Beam theory gives a free bar the fundamental \(f_1 = \dfrac{4.730^2}{2 \pi L^2} \sqrt{\dfrac{E I}{\rho A}}\) and its first overtones at 2.7565 and 5.4039 times it. The model lies close to it. Its overtones lie a little lower, because the bar also shears, which beam theory leaves out and which matters more for the shorter waves of the higher modes:

f, vectors, mesh, V = modes(bar(0))
f_beam = 4.7300407**2 / (2 * math.pi * L**2) * math.sqrt(E * t**2 / 12 / rho)
print("f1", f[0], "| beam", f_beam)
print("f2 / f1", f[1] / f[0], "| beam 2.7565 | f3 / f1", f[2] / f[0], "| beam 5.4039")

f1 552.6810885693504 | beam 555.4355406489958
f2 / f1 2.7315312224928654 | beam 2.7565 | f3 / f1 5.285645509802257 | beam 5.4039

The nodes of the fundamental are where its vertical motion along the middle of the bar changes sign:

mode = ngs.GridFunction(V)
mode.vec.FV().NumPy()[:] = vectors[:, 0]
xs = np.linspace(0, L, 3501)
up = np.array([mode(mesh(x, t / 2))[1] for x in xs])
print("nodes", xs[:-1][np.sign(up[:-1]) != np.sign(up[1:])], "| 0.2242 L", 0.2242 * L, 0.7758 * L)

nodes [ 78.4 271.5] | 0.2242 L 78.47 271.53000000000003

Tuning the overtone

A deeper arch lowers both frequencies, the fundamental more, so the ratio of the two grows with the depth. brentq finds the depth at which it is four (Figure 1):

def ratio(depth):
    f = modes(bar(depth))[0]
    return f[1] / f[0]


depths = np.arange(0, 9.5, 0.5)
ratios = [ratio(d) for d in depths]
depth = brentq(lambda d: ratio(d) - 4, 2, 9, xtol=1e-4)
f, vectors, mesh, V = modes(bar(depth))
print("depth", depth, "| f", f[:3], "| f2 / f1", f[1] / f[0], "| f3 / f1", f[2] / f[0])
# plt.plot(depths, ratios); plt.show()
# px.Plot(bar(depth))

depth 7.728074308615789 | f [ 225.05581979  900.2245988  2091.30820374] | f2 / f1 4.000005863628423 | f3 / f1 9.292397795770036

Figure 1: The overtone ratio against the depth of the arch

The arch takes 7.7 of the 13 in the middle and lowers the fundamental to about two fifths of the plain bar’s. The next overtone lands at about 9.3 times the fundamental, short of the ten that marimba makers aim for with a more shaped cut. Figure 2 shows the tuned bar and its first two modes.

(a) The tuned bar

(b) Fundamental, 225 Hz

(c) First overtone, 900 Hz

Figure 2: The tuned bar and its first two modes, the motion magnified

See netgen and NGSolve for meshing a planex2d face with netgen.