Measuring Absorbed Energy in the Human Auditory System Using Finite Element Models: A Comparison with Experimental Results

1. Introduction

The human auditory system (HAS) is an essential tool in performing numerous daily tasks. Therefore, it is important to be able to assess the proper functioning of the HAS in the most efficient way possible. The HAS is more complex than one might think at first glance, and there are still mechanisms whose exact functioning is unknown, or for which there are only hypotheses, such as the process of feedback from the inner and outer hair cells in the Basilar Membrane. In addition to this, there is the difficulty of accessing the HAS, which is embedded in the skull. Any intervention beyond the eardrum requires perforation or opening of the skull, with the corresponding consequences. These special circumstances have led medicine and engineering to develop non-invasive methods to determine the correct functioning of the HAS and, in case of malfunction, detect any associated pathology. The first and simplest test is audiometry, which involves subjecting the patient to the identification of beeps at different frequencies, increasing the intensity until surpassing the hearing threshold, at which point the patient informs the technician of the beep. Audiometry determines our hearing capacity in terms of frequency and decibels. Despite its simplicity, this test requires the patient’s cooperation and is not applicable to babies or individuals who are not conscious at the time. The measurement of impedance and/or Absorbed or Reflected Energy has been presented as an optimal solution.

Moller [1] first introduced the use of calibrated sound sources as a means for measuring impedances. The technique involves introducing a transmitter and receiver of waves into the auditory canal, so that the difference in energy between the emission and reception is evaluated. This difference can be related to the impedances of the different subsystems of the middle and inner ear, allowing for the inference of possible pathologies. This method, with some modifications or performances, is now extensively used. The following research [2,3,4,5] establishes an identical foundation in terms of the pressure source and receiver in the auditory canal. The difference between them lies in the technique used to approximate the calculation of eardrum impedance. Some of the studies from this time [6,7] begin to relate variations in impedance and, consequently, absorbed and reflected energies, to potential hearing loss.

In recent years, numerous studies have been conducted to improve the measurement of acoustic impedances [8,9,10,11,12], introducing Wideband Tympanometry (WBT) to replace the original Tympanometry. Nowadays, WBT is widely accepted and used for the calculation of impedances and reflected and absorbed energies. To do this, a passive test is performed that does not require any action from the patient. Some of these studies are more descriptive of the process and its proper setup [8,9], while others focus more on the relationship between the results and potential pathologies and hearing loss [10,11,12]. Norton or Thevenin equivalent sources are used for measuring device calibration, with only a single pressure measurement then being needed to calculate the external auditory canal (EAC) entrance impedance [13]. Other research has used the transfer function method, which relies on measuring the influence of the termination impedance on duct pressure [14], thereby requiring robust coupling of the device and the EAC [4,6]. Lanoye et al. [15] proposed a third method using impedance probes containing separated pressure and volume velocity sensors, although high levels of disturbance of the sound field at the measuring position are a problem when using this technique.

Most studies published to date have employed the two first methods where a sound (pressure) source and measuring device are usually placed inside the EAC. This is to avoid discontinuity between the device and the EAC to avoid exciting higher-order modes [16,17,18,19]. The frequencies of the minimum and maximum input impedances are mainly affected by the EAC length and its cross-sectional area. This means that the sound source and measurement device directly affect the impedance calculations because they change the canal length. One possible solution to this problem could be the application of inverse procedures to derive the EAC shape from its input impedance [12]. Thus, knowing the calculated cross-sectional area of the EAC allows the required eardrum impedance transformation for the energy calculations to be estimated.

To obtain middle ear diagnostics, the input impedances of the EAC measurements are transformed to represent the position of the eardrum and thus, record reflected energy (RE) and the energy absorbed, abbreviated as EA [5,7,8,13,14,15,16,17,18]. However, these impedance and energy results are often inaccurate. Thus, the main objective of this current work was to clarify the origin of possible errors in these measurements and their transformation as well as the consequences of these errors in terms of the impedances and energy calculated from them. In this sense, the use of finite element models (FEMs) is now widely accepted as an adequate and complementary research tool. Thus, this paper was based on several numerical simulations conducted using a FEM previously described and validated in the academic literature. This original work was first conducted for the outer and middle ear [20,21], was then applied to an inner ear FEM by employing a semiautomatic algorithm [22] and was finally coupled to the previous outer and middle ear FEMs [23]. Therefore, this FEM presents a complete fluid–structure interaction between the EAC, tympanic membrane (TM), and the oval window.

TM modelling was a crucial innovation with respect to previous models because the elements it used were better formulated and thereby eliminated problems associated with ‘shear locking’ elements used in thin membranes. Thus, together with adequate mesh convergence analysis, TM modelling provides sufficient guarantees of correct results. This means that, apart from geometric uncertainties caused by natural variability, most inaccuracy in this type of modelling comes from the difficulty of discerning the mechanical properties of some components such as the TM, tensors, and joints, among others. In this current article, we used assumed values for these components without attempting to discuss their accuracy. We simulated different combinations to discern the impact of each subsystem in the human auditory system (AS), building on a previous paper with a similar methodology [24] that determined how the AS influences pressure distribution in the EAC.

In this context, over the past century, Rong Gan has led a research group that has published several FEMs that have become an important source of inspiration for other researchers developing FEMs [25,26,27,28,29]. Gan’s work has focused on calculating EAC and eardrum impedances as well as the EA and ER from the TM [29]. Our work differs from that of Gan et al. in two clear ways. First, their main objective was to calculate the relationship between three middle ear disorders (otitis media, otosclerosis, and ossicular chain disarticulation), as simulated in their FEM, and any changes in the EA. Second, in contrast, our objective was to determine how the AS subsystems affect EA. We employed a range of FEMs, with the most basic one comprising an EAC and eardrum, to calculate the impedances and EAs. This methodology has already been applied successfully in previous work which gave us a better understanding of the mutual influence of each part of the AS on these factors.

It is important to understand the differences in our strategic goals and those of Gan et al. because these differing objectives affect the simulation conditions employed. Gan and colleagues aimed to offer a useful tool for the diagnosis of middle ear disorders by testing EAs through pneumatic otoscopic and wideband absorbance audiometry tests. Therefore, they literally reproduced the experimental setup boundary conditions and placed the sound source 20 mm from the eardrum to simulate the experimental test conditions in the FEM (see Figure 1). They determined that this approach gave them the best correlation between their numerical and experimental results. However, this does not mean that these numerical results most closely resemble reality. Indeed, controversy remains regarding the accuracy of the experimental methods used to measure impedance and in turn, the subsequent energy calculation results [30,31,32,33,34]. Therefore, more theoretical background work will still be required to increase our knowledge in this area.

2. Theoretical Background

The EA is calculated based on the ER. It describes the fraction of incident acoustic power reflected by the TM, where a reflective power of 1 corresponds to complete reflection of all the acoustic power and a reflectance of 0 corresponds to the condition in which all the power is absorbed by the TM [35]. First, the characteristic impedance of the EAC was calculated as:

$$Z_c=\frac{\rho c}{S}$$

(1)

where $\rho$ is the density of the air contained in the EAC, c is the speed of sound in air, and S is the cross-sectional area of the EAC. In this study, air density and sound speed values of 1.21 kg/m³ and 343 m/s were used, respectively. To calculate the TM impedance, the impedance of each of the elements of the FEM for the TM was determined as a function of the acoustic pressure $p_{TM}^i\left(f\right)$, with $v_{TM}^i\left(f\right)$being the velocity in element i of the TM, and where $A_{TM}^i$ was the element area. Both the velocity and pressure depend on the frequency f, so we used values obtained for a frequency range of 0.1–10 kHz according to the FEM:

$$Z_{TM}^i\left(f\right)=\frac{p_{TM}^i\left(f\right)}{A_{TM}^i{v}_{TM}^i\left(f\right)}$$

(2)

The total TM impedance, $Z_{TM}\left(f\right)$, was obtained by adding the impedances of each of the elements of the TM surface in parallel:

$$Z_{TM}\left(f\right)=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\sum _i\left[\frac{1}{Z_{TM}^i\left(f\right)}\right]$}\right.$$

(3)

The acoustic impedance in the EAC was calculated based on $Z_{TM}\left(f\right)$ and $Z_c$:

$$Z_{EC}\left(f\right)=Z_c\frac{Z_{TM}\left(f\right)+j{Z}_c\tan \left(kL\right)}{Z_c+j{Z}_{TM}\left(f\right)\tan \left(kL\right)}$$

(4)

where k is the wave number and L is the distance between the TM and the location of the measurement points in the EAC, which was at 30 mm in our study. The reflected acoustic pressure is obtained using the expression:

$$R_P\left(f\right)=\frac{Z_{EC}-Z_C}{Z_{EC}+Z_C}$$

(5)

Thus, the ER is calculated based on the reflected acoustic pressure $R_P\left(f\right):$

$$ER\left(f\right)={\left|R_P\left(f\right)\right|}^2={\left|\frac{Z_{EC}-Z_C}{Z_{EC}+Z_C}\right|}^2$$

(6)

Finally, the EA was obtained as a function of the frequency calculated as:

$$EA\left(f\right)=1-ER\left(f\right)$$

(7)

3. Materials Models

Finite Element Model

The construction and validation of the FEM applied in this section, as well as its material properties, have been previously published [20,21,22,23,24]. However, in this current work, the EAC was re-meshed with hexahedral elements [36,37] to help improve the accuracy of the FEM post-process calculations (Figure 2).

As summarised in Table 1, three different FEM combinations were simulated to discern the impact of the different AS subsystems.

4. Results

The TM was modelled using 7880 hexahedral elements, with each one itself comprising 8 nodes, although only the 4 that formed the surface in contact with EAC air were valid in our calculations. Thus, a total of 8193 nodes formed the TM surface. With harmonic analysis, we obtained 16,386 real and imaginary values each for pressure and speed. Thus, when we performed the harmonic calculation for the 0.1–10 kHz range, we produced matrices comprising a total of 1,638,600 data points each for pressure and speed. In turn, the area was composed of only the 7880 element data points. We then used ANSYS 2022R2 engineering simulation software to obtain a total of 3,285,080 values for use in impedance and EA calculations using MATLAB R2023A software.

4.1. Tympanic membrane velocity

The average of all the elements in the TM velocity module is shown in Figure 3. There were three clearly differentiated resonance zones. One was at around 800–1000 Hz because of the eardrum itself, while the other two, located at 4 kHz and 9 kHz, were the result of the characteristic resonances of the EAC [24]. These results highlight the importance of modelling the ossicular chain, as there is a significant difference in the outcomes when the ossicular chain and cochlea are not modelled. On the other hand, there are no significant differences observed between modelling the simplified cochlea and the realistic spiral model, resulting in significant computational savings.

4.2. External auditory canal characteristic impedance

The characteristic impedance of the EAC, Z_c depends on the air density ($\rho$), speed of sound in air (c), and cross-sectional area of the canal (S). The cross-sectional area was measured at the entrance of the canal where the pressure source was located and was 89,652 mm². The density and speed of sound values were 1.21 kg/m³ and 343 m/s, respectively, and so the characteristic impedance of the EAC was 4.6294 × 10³ Pa s/mm³. This value was used for all the FEMs and was frequency independent.

4.3. Tympanic membrane impedance (direct calculation)

Figure 4 represents the module and phase TM impedance in a frequency range of 0.1 to 10 kHz. Figure 4A shows a comparison of the module impedance for the three studied FEMs and showed a decrease in the impedance with frequency up to 800–1000 Hz, after which it increased. This minimum coincided with the first eardrum resonance frequency in that range. Figure 4B shows how the system formed by the EAC and TM opposes the system of least resistance to wave propagation and so was where the highest speeds occurred. The phase impedance started at a value of −80° for 100 Hz and increased with the frequency until it reached 0° for a range of 700–900 Hz. For frequency values between 8000–9000 Hz, the phase impedance reached a maximum between 60–80°. Therefore, the phase advanced until it reached the resonance frequency at around 900 Hz, with the phase also changing from lagging to leading at this point. In addition, there was a tendency towards an asymptotic phase shift at 180°, meaning that this was a second order system.

4.4. External auditory canal impedance (backward calculation) based on tympanic membrane impedance.

Figure 5A shows the module EAC impedance for each of the studied FEMs, which strongly coincided, especially at high frequencies. The EAC curves presented two minima, one at 3000 Hz and the other at 9000 Hz, with a maximum at around 6000 Hz. In turn, Figure 5B shows the phase EAC impedance curves for each FEM, which also strongly coincided at high frequencies. The TM stopped influencing the results in both the module and phase EAC impedance at around 1000 Hz [40,41], with curve disturbances at a lower range than those present at 3000, 6000, or 9000 Hz. The first and second resonance frequencies of the EAC were at 3000 and 9000 Hz, respectively, while the anti-resonance of the canal was at 6000 Hz.

4.5. Energy absorbance at the tympanic membrane

The results of EA for our models are presented In Figure 6, the EA curve presented a maximum for a frequency value around 700–900 Hz in all three tested FEMs, reaching a minimum for low frequency values between 100–200 Hz and a maximum at 9000–10,000 Hz. Figure 6 shows how the maximum EA values coincided with the minimum of the TM impedance, as shown in Figure 4B.

In Figure 7 are shown the EA results published. Zhang [29] presented results from numerical simulations by FEM. Other results [5,35,43,44] belongs to experimental tests. These experimental tests were carried out placing the pressure source inside the EAC, changing the natural EAC boundary conditions. There are differences between the experimental results, but It could be separate in two groups: on the one hand, Feeney [45] has the maximum before 1 kHz and fast decay; on the other hand, rest of experimental results have the maximum at 4 kHz. Our FEM results are more coincident with those presented by Feeney [45]

5. Discussion

To calculate the impedances, it appears that a well-defined cochlea model is unnecessary, especially at high frequencies. To validate this finding, we compared our data with those provided in previous academic publications. As shown in Figure 8 for module (A) and phase (B) EAC impedance, the values obtained with our FEMs corresponded to the data obtained by Zhang and Gan [36]. The blue curve represents the impedance of the EAC calculated according to Zhang and Gan to which a correction method [42] was applied to simulate our test condition (L = 0.03 m) and our EAC section (S), which thereby modified the value of the characteristic impedance (Z_c). These results strongly concurred with the complete FEM model (CATIOSCO) for this correction.

6. Conclusions

Here we established a process for the numerical calculation of impedances and acoustic EA by the TM. This process starts with post-processing in ANSYS software followed by data exportation and processing in MATLAB. Regarding the impedance of the EAC, it is worth noting that the results obtained in this study agreed with those published by Gan and colleagues once the length of the canal was modified to 0.03 m, with values ranging from 10⁻⁴ to 10⁻¹ Pa·s/mm³, valleys at 3 and 9 kHz, and a peak at 7 kHz. Modelling of the ossicular chain and cochlea were shown to be crucial to determine the impedance of the TM. Without these elements, the measured impedance dropped to 0.6·10⁻² Pa·s/mm³ with a minimum around 800 Hz, while modelling the complete system reduced the minimum to 2·10⁻² Pa·s/mm³, shifting the valley frequency to approximately 1 kHz. Of note, the experimental values previously reported [2,5,36] closely matched those obtained in this current study.

In this study were significantly different when the ossicular chain and cochlea were introduce ed into the models. Without these elements, EA reached a value of almost 1 at around 900 Hz. This is quite logical because this frequency coincides with the first natural frequency of the eardrum at which almost all incoming energy is absorbed. When the ossicular chain and cochlea were modelled, this maximum value reduced to 0.6 and the frequency slightly increased to 1 kHz. Nonetheless, experimental results showed discrepancies in these figures. Previous work [46] exhibited a similar maximum EA frequency but with a value close to 1 while other results [5,23,25,26,36] showed that the first part was nearly identical to our findings in this current study, with a maximum around 1 kHz and a value of 0.6, except that the EA continued to increase to 3–4 kHz, with values between 0.7 and 1.

Considering the equation used to calculate EA (Equation 7) and referring to its root equations (Equations 4– 6) involving Z_TM, Z_C, and Z_EC (the TM and EAC impedances), we first showed that the Z_TM we obtained coincided with the findings from experimental work presented elsewhere. Thus, we deduced that Z_EC caused the differences between the modelled and experimental results mentioned above. This was because placement of the measuring device in the EAC considerably reduced its length, thereby affecting these impedance values, as demonstrated in Figure 8. This may be because it is difficult to establish the conditions of a ‘normal ear’ as a FEM [23].