Acoustic simulation with 1D Finite Element Method Analysis¶

This algorithm computes the impedance spectrum of a didgeridoo bore using 1D Finite Element Analysis. This simulation method should be more precise than the Transmission Line Model, but it has a higher computational complexity. The artificial evolution computes 1000s and 10000s of these simulations over multiple hours, so I prefer the faster Transmission Line Model. Below is the source code and a brief explanation of the mathematics.

In [7]:

import numpy as np
import skfem as fem
from skfem.helpers import dot, grad
import matplotlib.pyplot as plt

c = 343000.0  # Speed of sound in mm/s

def fem1d(geo, frequencies):


    x_coords = geo[:, 0]
    diameters = geo[:, 1]

    def get_area(x):
        d = np.interp(x, x_coords, diameters)
        return np.pi * (d / 2.0)**2

    # 2. Mesh and Basis
    mesh = fem.MeshLine(np.linspace(0, x_coords[-1], 600))
    element = fem.ElementLineP1()
    basis = fem.Basis(mesh, element)

    # 3. Define Forms
    def stiffness_fun(u, v, w):
        return get_area(w.x[0]) * dot(grad(u), grad(v))

    def mass_fun(u, v, w):
        return get_area(w.x[0]) * u * v

    K = fem.BilinearForm(stiffness_fun).assemble(basis)
    M = fem.BilinearForm(mass_fun).assemble(basis)

    # 4. Boundary Setup
    # Bell end: Pressure = 0 (Dirichlet)
    bell_dof_data = basis.get_dofs(lambda x: np.isclose(x[0], x_coords[-1])).nodal
    # Extract indices from dictionary if necessary
    bell_dofs = np.concatenate([v for v in bell_dof_data.values()]) if isinstance(bell_dof_data, dict) else bell_dof_data

    # Mouthpiece: Setup the source vector
    b_mouth = np.zeros(basis.N)

    # FIXED: Safely extract nodal indices from the Dofs object (handling dict structure)
    mouth_dof_data = basis.get_dofs(lambda x: np.isclose(x[0], 0)).nodal
    if isinstance(mouth_dof_data, dict):
        mouth_indices = np.concatenate([v for v in mouth_dof_data.values()]).astype(int)
    else:
        mouth_indices = np.array(mouth_dof_data).astype(int)

    b_mouth[mouth_indices] = 1.0 

    # 5. Frequency Sweep
    impedance_magnitude = []

    # Add viscothermal damping to prevent infinite peaks and simulate wall friction
    # damping factor alpha is proportional to sqrt(frequency)
    for f in frequencies:
        omega = 2 * np.pi * f
        k = (omega / c) - 1j * (2e-6 * np.sqrt(f)) 
        k_sq = k**2
        
        A = K - k_sq * M
        
        # Solve for complex pressure p
        p = fem.solve(*fem.condense(A, b_mouth, D=bell_dofs))
        
        # Input Impedance Z = p_mouth / U_mouth. Since U=1, Z = p[mouth]
        z_val = p[mouth_indices[0]]
        impedance_magnitude.append(np.abs(z_val))

    impedance_magnitude = np.array(impedance_magnitude)
    return impedance_magnitude

geo = np.array([[0, 30], [1000, 30], [1400, 64]])
freq_grid = np.linspace(30, 1000, 1000)
z = fem1d(geo, freq_grid)
z = np.log2(z)

plt.plot(freq_grid, z)
plt.xlabel("Frequency")
plt.ylabel("Impedance")

Out[7]:

Text(0, 0.5, 'Impedance')

No description has been provided for this image

1D Finite Element Method¶

1. The Governing Equation: Helmholtz Equation¶

The code models the didgeridoo using the Webster Horn Equation, which is a 1D version of the Helmholtz equation that accounts for a changing cross-sectional area $S(x)$. The physics follows this second-order differential equation for acoustic pressure $p$: $$\frac{1}{S(x)} \frac{d}{dx} \left( S(x) \frac{dp}{dx} \right) + k^2 p = 0$$ $S(x)$: Cross-sectional area at position $x$. $k$: The complex wave number ($k = \frac{\omega}{c}$). $p$: Acoustic pressure.

2. The Weak Form (Finite Element Formulation)¶

To solve this with FEM, the code converts the differential equation into a "weak form" by multiplying by a test function $v$ and integrating over the length of the tube $L$. After applying integration by parts, we get: $$\int_0^L S(x) \frac{dp}{dx} \frac{dv}{dx} dx - k^2 \int_0^L S(x) p v dx = 0$$ In the code, this is split into two bilinear forms: Stiffness Matrix ($K$): get_area(w.x[0]) * dot(grad(u), grad(v)) — Represents the spatial variation of pressure. Mass Matrix ($M$): get_area(w.x[0]) * u * v — Represents the "storage" of acoustic energy in the volume.

3. Discretization and Linear Algebra¶

The code uses skfem to divide the instrument into 600 small linear elements (ElementLineP1). This turns the continuous calculus problem into a discrete matrix equation: $$(K - k^2 M) \mathbf{p} = \mathbf{b}$$ $K$ and $M$ are global matrices assembled from the segment geometries. $\mathbf{p}$ is the vector of unknown pressures at each node. $\mathbf{b}$ is the "source" vector (the input from your lips at the mouthpiece).

4. Complex Wave Number and Damping¶

In an ideal world, $k = \omega / c$. However, without damping, the resonances would be infinitely sharp (mathematically singular). The code introduces an imaginary component to $k$ to simulate viscothermal losses: k = (omega / c) - 1j * (2e-6 * np.sqrt(f)) This imaginary term represents energy absorbed by the walls of the didgeridoo. The $\sqrt{f}$ dependency is mathematically accurate for boundary layer losses in acoustic tubes.

5. Boundary Conditions¶

The code defines how the sound behaves at both ends: The Mouthpiece ($x=0$): It sets a source term b_mouth[mouth_indices] = 1.0. This effectively simulates a unit volume velocity input. The Bell ($x=L$): It sets p = 0 (Dirichlet condition) using bell_dofs. In physics, this is an "open end" approximation where pressure drops to zero because the tube meets the infinite atmosphere.

6. Calculating Impedance¶

Acoustic Impedance $Z$ is defined as the ratio of Pressure $p$ to Volume Velocity $U$: $$Z = \frac{p}{U}$$ Because the code sets the input source ($U$) to a constant $1.0$, the resulting pressure value at the first node (p[mouth_indices[0]]) is the impedance. The impedance_magnitude is simply the absolute value of this complex result.