In the Ear Canal - Part 1 episode we built a first ear canal model. We used a freely available cast of an ear canal and manipulated its geometry. We obtained some initial results, but it was clear that the accuracy left much to be desired. In this episode we step back: we solve an eigenproblem on a better geometry and compare the results with a benchmark problem.

Project Files

All the files used for this project are available at the repositories below:

• Ear Canal 2.

The Problems of the Previous Model

The previous model had multiple issues, including geometry, governing equation and post-processing techniques. From a geometry point of view, it used an ear canal cast from a particular person. While this is not wrong, the specific of this person anatomy might make comparison with data from literature harder, as that typically involves averages over multiple people. This episode focus on a better geometry.

A Better Geometry

It is not straightforward to say what a better geometry is. Ear canals, like most anatomical parts, are subject to large variation. A good geometry definition could be that of a reasonable average over a large number of diverse humans. However, how to measure the ear canal anatomy, and then reasonably average the results, is a complex problem on its own.

Thankfully we do not have to do that step. The ITU-T P.57 recommendations contain specifications for various types of ear simulators, including one that is anatomically shaped: the Type 4.3. It is the specifications of this simulator that we will use to define our geometry.

Of course, the ITU methodology (which you can read in the standard) is not immune to critique. Brian Johansen, Per Rasmussen, and Morten H. Pedersen have made a veiled one in an article they released for audioxpress:

In addition to the attempt to create a “one size fits all” ear simulator that, in reality, fits no one but simply represents an “average ear” geometry, there is a rarely explored dimension at play.

However, using an ITU standardized geometry has the following pro: our results can be compared with other people doing similar work. For example:

• The Digital Twin of a New and Standardized Fullband Ear Simulator by Lars Birger Nielsen1 and Mads Herring Jensen
• Type 4.3 Ear Simulator from the COMSOL Application Gallery.

It is therefore a better starting point.

Implementing ITU-T P.57

Difficulties

The implementation of ITU-T P.57 was much harder than I anticipated. If you read ITU-T P.57 (06/2021) you will see that Table 6 specifies an ear canal center line. Along this center line various cross sections of the ear canal are presented in planes orthogonal to the center line. The unit vectors defining these planes are presented in Table 7, while the cross sections shapes are presented in Figure 12.

At first sight is all seems simple: sample the cross sections outlines and put the points in space. Fit some closed loop curve through the points of each section outline and then loft a solid through the loops.

This operation can be in fact done with Gmsh by:

• Defining Points along the outline of each cross section.
• Define a closed BSpline through them, and add each one as a Curve loop.
• Create a Thrusections between the loops.

However, despite my best efforts, this very often results in a geometry where only an outer surface is present, without a solid. This happens, I think, due to the outer surface curling excessively in certain spots. And this is, in turn, created by the step by which the cross sections are provided, together with some sharp turns of the ear canal center line.

My attempts at creating a good ear canal geometry are perhaps too lengthy to discuss in here. If you are interested, you can read the code at acoupy_ears. This is the python package that I written to implement the geometry and it will be better documented in the future. As of now, I could only produce:

• A flat terminated ear canal model ending at 28 mm (Ear Canal End).
• A flat terminated ear canal model ending at the last EEP (Ear canal Entrance Point) projection point at 32.5 mm (EEP End).

If you want to install acoupy_ears into any python environment you can use:

pip install git+https://gitlab.com/acoupy/acoupy_ears.git

After installation you will have the command build available. Use build --help to see the help. The sub-commands allow you to save a BREP file of the geometry on disk. You can process this file with Salome as per usual.

The Resulting Geometry

The ear canal geometry up to its end at 28 mm is shown below. Some of the unnatural bends of the surface are visible. Thankfully, a mesh size that works up to 20 kHz is still about 1 mm. This means that a mashing algorithm will still be successful in meshing this geometry, since the nodes end up falling on opposite sides of the unnatural bends.

The same holds for the ear canal geometry up to the last EEP projection point at 32.5 mm, which is shown below (on the same scale of the figure above).

Both geometries can be loaded into Salome and exploded in a solid. This solid can then be further exploded in 3 faces: the two flat terminations and the “tube” surface. The solid can therefore be meshed as usual with convenient mesh groups for boundary conditions. If this doesn’t make any sense you can refer to a previous tutorial, for example Home Studio - Part 1 | Meshing, to see how you can load a BREP file into Salome, process its geometry and mesh it.

In this episode both ear canal models have been meshed with NETGEN to a maximum (and minimum) size of 1 mm. For details you can refer to the Salome study file.

The Benchmark Problem

I decided to introduce a simple benchmark to evaluate our ear canal models accuracy. As a first approximation, an ear canal is a cylinder. However, various sources differ on the size of this cylinder. This confusion is perhaps due to the complexity in measuring the length of an ear canal: where does exactly the ear canal start?

For example, Psychoacoustics Facts and Models by Hugo Fastl and Eberhard Zwicker states (Third Edition, page 23):

It [the ear canal] acts like an open pipe with a length of about 2 cm corresponding to a quarter of the wavelength of frequencies near 4 kHz.

In Fundamentals of Acoustics by Lawrence E. Kinsler, Austin R. Frey, Alan B. Coppens and James V. Sanders instead we read (Fourth Edition, page 312):

The auditory canal is an approximately straight tube, about 0.8 cm in diameter and 2.8 cm long, closed at its inner end by the tympanic membrane (eardrum). The lowest resonance of this tube is broadly peaked around 3 kHz and affords appreciable gains from about 2 kHz to 6 kHz.

Wikipedia, citing Structural and functional anatomy of the outer and middle ear by Faddis, B. T. (which I couldn’t access) states:

The adult human ear canal extends from the auricle to the eardrum and is about 2.5 centimetres (1 in) in length and 0.7 centimetres (0.3 in) in diameter.

Clearly, there are many different values floating around, both for the geometry and the fundamental resonance frequency. Moreover, the fundamental resonance frequency is not only affected by the geometry, but also by the boundary condition. There is a difference between considering, for example, the ear canal geometry as fully enclosed by a rigid boundary versus terminated at one side by the tympanic membrane impedance and open on the other. The sources above do not specify exactly in what condition the ear canal resonance has been measured or modeled (open at its input or driven at its input?). We will see that we can piece together something, however.

In this article we will do the simplest possible model. We will follow Kinsler and Frey and use 0.8 cm for diameter and 2.8 cm for height. This because the 2.8 cm length reported by them ties nicely with length of one of our ITU-T P.57 models. And we will solve for the resonances of the ear canal models with a rigid boundary condition on every boundary. That is, we consider the normal modes of the air contained in the ear canal when this is plugged by a rigid flat cap, assuming that the inside of the ear canal, including the ear drum, is a rigid boundary too.

We can compute the normal modes of the cylinder analytically in this case. Therefore, we will be able to assess whether Elmer gives us good answers and, if so, we will switch to the ITU-T P.57 geometries with confidence.

The benchmark mesh, also meshed to 1 mm size, is shown below in comparison with the ear canal meshes from ITU-T P.57.

Note that in the figure above the benchmark cylinder is not in the same geometrical scale as the other models as its geometry is defined in a different reference frame.

The Models

In this episode I will discuss the results from Elmer’s Wave Equation eigensolver for all the discussed geometries. That of setting up this type of simulation with Elmer is a well discussed topic on this website. If you are not familiar with it you can refer to the Rigid Walled Room Series. In here I will simply discuss how the problem has been setup in broad strokes. This type of problem is actually easier to setup from scratch than through ElmerGUI. Below is the summary of steps to setup the problem.

Once you have done meshing your geometry in Salome you can export the mesh as a UNV file (see, for example, Rigid Walled Room | Meshing). Once that is done, you can use ElmerGrid to convert the UNV mesh file into a mesh in Elmer format:

ElmerGrid 8 2 my_mesh.unv -autoclean

This will create a new my_mesh folder with the various files defining the mesh according to Elmer’s needs. Now you can create SIF file for the problem. You can call you SIF file case.sif. If you put it in the same folder as the newly create my_mesh folder you will need it to have this Header:

Header
  CHECK KEYWORDS Warn
  Mesh DB "." "my_mesh"
  Include Path ""
  Results Directory ""
End

Note how the Mesh DB field points to the mesh folder.

This is a steady state simulation and we must not forget about Coordinate Scaling, if needed. As we learned in the Rigid Walled Room Series it is also very useful to have Elmer outputting a series of VTU files, one for each eigenvalue-eigenfunction pair. So, the Simulation section would look like this:

Simulation
  Max Output Level = 5
  Coordinate System = Cartesian
  Coordinate Mapping(3) = 1 2 3
  Simulation Type = Steady state
  Steady State Max Iterations = 1
  Output Intervals(1) = 1
  Coordinate Scaling = 0.001
  Solver Input File = case.sif
  Post File = case.vtu
  vtu: Eigen Analysis = Logical True
End

The solvers and equations sections are the most important. As is tradition for these problems we can use two solvers: one to actually solve the problem and one to save the eigenvalues to file. For the former, we need to specify we are doing an solving an eigenvalue problem and provide the relevant analysis parameters. For higher accuracy we should also use second order P-elements. This is shown below:

Solver 1

  ! Equation
  Equation = "Wave Equation"
  Procedure = "WaveSolver" "WaveSolver"
  Variable = "Pressure"
  Variable DOFs = 1
  
  ! Solver
  
  ! Steady State
  Steady State Convergence Tolerance = 1e-09
 
  ! Linear System
  ! Linear System Scaling = Logical True
  ! Linear System Symmetric = Logical True
  Linear System Solver = string "Iterative"
  Linear System Direct Method = Umfpack
  Linear System Iterative Method = BiCGStabl
  Linear System Max Iterations = integer 10000
  Linear System Convergence Tolerance = real 1e-8
  BiCGstabl polynomial degree = 4
  Linear System Preconditioning = ILUT
  Linear System ILUT Tolerance = 1.0e-3
  Linear System Abort Not Converged = True
  Linear System Residual Output = 50
  Linear System Precondition Recompute = 1
  
  ! Special
  Eigen Analysis = True
  Eigen System Damped = False
  Eigen System Values = 10
  Eigen System Select = smallest magnitude
  Eigen System Convergence Tolerance = Real 1.0e-9
  Eigen System Max Iterations = 10000
  ! Eigen System Shift = Real 1000
  Eigen System Normalize To Unity = Logical True
  Element = "p:2"
  
End

Solver 2
  Equation = "SaveScalars"
  Filename = eigenvalues.dat

  Procedure = "SaveData" "SaveScalars"
  Save EigenValues = True

  ! Compare the norm of the 10th eigenvalue
  Show Norm Index = 10 
  
End

Equation 1
  Name = "Wave Equation"
  Active Solvers(1) = 1
End

The rest of the SIF setup is pretty straightforward. You will need to define a Material and a Body so that both the material and Equation 1 are assigned to it. For the materials you need to ensure the parameters needed by the solver are provided (Sound Speed, Sound Damping and Sound Reaction Damping).

The boundary conditions are not essential for eigenvalue problems: Elmer will apply a rigid one pretty much no matter how you set them. Some insights on why are given in this forum post. So, not only a fully rigid boundary is the simplest problem that we can treat. When it comes to eigenvalue problems Elmer’s will force our hand in this regard.

For SIF files examples you can refer to the repository.

To finish up the setup create the ELMERSOLVER_STARTINFO with case.sif as its contents. Now you can run ElmerSolver from the SIF file folder and have Elmer running your simulation.

The Solution

The Cylinder

The cylinder is our benchmark problem, and for this we have closed form solutions. Fundamentals of Acoustics treats the problem of a rigidly enclosed cylindrical cavity succinctly (Fourth Edition, page 249). From it we have:

$$ f_{l,m,n} = \frac{c}{2 \pi} \sqrt[+]{ \left( \frac{j_{m,n}^{\prime}}{a} \right)^2 + \left( \frac{l \pi}{L} \right)^2 } $$

Where $c$ is the speed of sound (a material parameter), $a$ is the radius of the circular base of the cylinder, $L$ is the height of the cylinder and the $j_{m,n}^{\prime}$ values are the $n$-th extrema of the $m$-th Bessel function of the first kind. $l,m,n$ are 3 positive integers (zero included) that parameterize the mode. As for the rigid walled room, the cylinder is a 3D system and its normal modes are sustained by reflections from the boundaries. Therefore, these integers are the numbers of null surfaces along each of the reference system directions. However, for the cylinder we (of course) use a cylindrical coordinate system. This means that $l,m,n$ are the number of nulls in the axial, angular and radial directions respectively. The $j_{m,n}^{\prime}$ values are the most complex value to obtain for this equation. There must be plenty of way to get those by using any technical computing language, like Python or Julia. But also, they can be found tabulated in literature. For this article I will be using Extrema for Bessel Functions of the First Kind, J.

As a tip, if you like to process data with Python, the TSV file that you can download from the link above can be loaded quite nicely as a Polars dataframe once you extract the archive:

import polars as pl

j_df = pl.read_csv('BesselMaxAndMin.tsv', has_header=False, separator='\t')

As for the Elmer solution, remember that Elmer computes the eigenvalues $\lambda_k$, with $k$ an integer index counting among the eigenvalues that Elmer found. The associated modal frequency is given by:

$$ f_k = \frac{\sqrt[+]{\lambda_k}}{2 \pi} $$

How the $\lambda_k$ (and therefore the $f_k$) are sorted depends on the solver settings. In this case we are asking Elmer for the 10 smallest eigenvalues, sorted in ascending order. Therefore, we can go through the sorted $f_{l,m,n}$ values with the same index $k$.

The table below compares the first modal frequencies with frequency lower than 20 kHz (the higher frequency for our mesh).

The $k=0$ mode is, as always for the first mode, a “DC mode”. In other words, a constant pressure value. We can discard this mode (in fact, we could have done so by adding a small Eigen System Shift value in the solver settings). Overall, the agreement between Elmer and the exact solution is very good. According to Psychoacoustics Facts and Models (Third Edition, page 183) the minimum just noticeable frequency variation for humans is 3.6 Hz. The error between Elmer and the exact solution is lower than this for all the modal frequencies, showing that we can consider Elmer accurate for psychoacoustics purposes. Since we have found good agreement with the benchmark problem we can move on to the more complex geometry.

But first, a note. You have probably noticed that the first mode is 6125 Hz. This is quite a bit different from the values reported by various sources such as the ones listed above. This must be happening because the resonance values for the ear canal given by those sources are not for a fully enclosed ear canal. By taking a look at Fundamentals of Acoustics (Fourth Edition, page 272) we can see that different formulae hold for different boundary conditions at the ends of a cylindrical pipe. Most likely the resonances listed above were for the case in which a closed pipe is driven at one end by an uniform velocity source. In fact, in this case:

$$ f_n = \left(2n - 1 \right) \frac{c}{4 L} $$

With $n$ and integer index starting at 1. This would provide a first resonance at 3062.5 Hz, pretty consistent with the resonance reported above by the same source. This is a nice reminder of how critical boundary conditions are to the behavior of a system!

The Ear Canal Models

The table below compares the resonance frequencies of the two ear canal models. Clearly, the one with the most similar results to the benchmark problem is the 28 mm, the same length as the benchmark cylinder. The ear canal geometry has lower volume with respect the cylinder. It is therefore to be expected for the resonance frequencies to be all slightly higher with respect the benchmark. The longer system, the one extending to the EEP projection point, has a lower fundamental due to its larger volume. Therefore, it also has a higher number of resonance frequencies below 20 kHz.

We can now take a look at the modal patterns. Recall that the modal patterns are dimensionless, and normalized to take values between -1 and 1.

For the 28 mm model we can compare directly with the benchmark (note that in the figures the geometries are not on the same scales due to the different reference frames between the two models). In this case we can see that all the modal patterns are longitudinal, as they develop along the axis of the canal.

For the 32 mm model we have, instead, some more complex patterns.

Conclusion

In this episode we built a model of the ear canal based on ITU-T P.57. We also built a simple cylindrical ear canal benchmark model. We used these models to compute the resonance frequencies and modal patterns of the volume of air enclosed by the various geometries, with all rigid boundaries. We found that Elmer’s produces very good agreement of the benchmark model with theory. In addition, the results for the 28 mm ITU-T P.57 model are similar to that of the benchmark model, which makes sense given the matched lengths of these geometries.

However, we also found that the ear canal resonances are largely affected by boundary conditions, especially at the ear canal opening. In the following episode we will take these models and produce a frequency sweep for a uniform velocity source at the ear canal input to study resonances and frequency response when driving the canal and its opening.

License Information

This work is licensed under a Creative Commons Attribution 4.0 International License.

• Physics
• Elmer
• FEM
• Linear System
• Single Physics
• Series: Ear Canal