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

1. Introduction

In recent years, numerous studies have been developed to measure acoustic impedances[1,2,3,4,5,6,7,8,9,10,11]. Moller [12] first introduced the use of ‘calibrated sound sources’ as a means for measuring impedances, a method that is now extensively used [13,14,15,16,17,18,19,20]. 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[21]. Other research has used the transfer function method, which relies on measuring the influence of the termination impedance on duct pressure [22], thereby requiring robust coupling of the device and the EAC [16,17]. Finally, Lanoye et al. [23] 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 [24,25,26,27]. 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, 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 [35,36], was then applied to an inner ear FEM by employing a semiautomatic algorithm [37], and was finally coupled to the previous outer and middle ear FEMs [38]. 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 [39] 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 [40,41,42,43,44]. Gan’s work has focused on calculating EAC and eardrum impedances as well as the EA and ER from the TM [44]. 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. 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 [45,46,47,48,49,50]. 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 [31]. First, the characteristic impedance of the EAC was calculated as:

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

(1)

where ρ 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. Finite element model

The construction and validation of the FEM applied in this section, as well as its material properties, have been previously published [35,36,37,38]. However, in this current work, the EAC was re-meshed with hexahedral elements [51,52] 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 7,880 hexahedral elements, with each one itself comprising 8 nodes, although only the 4 that formed the surface in contact with CAE air were valid in our calculations. Thus, a total of 8,193 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 7,880 element data points. We then used ANSYS engineering simulation software to obtain a total of 3,285,080 values for use in impedance and EA calculations using MATLAB software.

A. 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–1,000 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 CAE [53].

B. External auditory canal characteristic impedance.

The characteristic impedance of the EAC Z_c depends on the air density, ρ, 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 we​are 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 FAMs and was frequency independent.

C. 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–1,000 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 CAE and TM opposes the system of least resistance to wave propagation and so was where the highest speeds occurred. Figure 4B shows how 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 8,000 -9,000 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.

D. 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 3,000 Hz and the other at 9,000 Hz, with a maximum at around 6,000 Hz. In turn, Figure 5B shows the phase EAC impedance curves for each FEM, which also strongly coincided at high frequencies. The WM stopped influencing the results in both the module and phase EAC impedance at around 1,000 Hz, with curve disturbances at a lower range than those present at 3,000, 6,000, or 9,000 Hz. The first and second resonance frequencies of the EAC were at 3,000 and 9,000 Hz, respectively, while the anti-resonance of the canal was at 6,000 Hz.

E. Energy absorbance at the tympanic membrane

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 9,000–10,000 Hz. Figure 6A shows how the maximum EA values coincided with the minimum of the MV impedance, as shown in Figure 4B.

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 [44]. The blue curve represents the impedance of the EAC calculated according to Zhang and Gan to which a correction method [50] 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 External Auditory Canal, it is worth noting the agreement of the results obtained in this study with those of GAN once the length of the canal is modified to 0.03 m, with values ranging from 10^-4 to 10^-1 Pa·s/mm³, with valleys at 3 and 9 kHz, and a peak at 7 kHz. The modelling of the ossicular chain, as well as the cochlea, are shown to be crucial in the impedance of the Tympanic Membrane. Without these elements, the impedance drops to 0.6·10^-2 Pa·s/mm³ with the minimum around 800 Hz, while modelling the complete system reduces the minimum to 2·10^-2 Pa·s/mm³, shifting the valley frequency to approximately 1 kHz. The experimental values [13,15,44] closely match those obtained in this study.

Regarding the values obtained for EA (Energy Absorption), there is a significant difference in this study when the ossicular chain and cochlea are introduced. Without these elements, EA reaches a value of almost 1 at around 900 Hz, which is quite logical as this frequency coincides with the first natural frequency of the eardrum, absorbing almost all incoming energy. When modelling the ossicular chain and cochlea, this maximum value is reduced to 0.6, and the frequency slightly increases to 1 kHz. Experimental results show discrepancies; This work [11] exhibits a similar frequency of maximum EA but with a value close to 1. Other results (15,31,33,34,44) have a first part nearly identical to that of this study, with a maximum around 1 kHz and a value of 0.6, but the significant difference is that EA continues to rise to 3-4 kHz with values between 0.7 and 1. Considering the equation used to calculate EA (Equation 7) and referring to the previous equations (Equations 4-6), it is deduced that the terms ZTM, ZC, and ZEC (Tympanic and Auditory Canal Impedances), once it has been shown that ZTM coincides between experimental work and the one presented in this study, it is deduced that ZEC is what causes the difference in results. This is due to the placement of the measuring device within the canal, considerably reducing its length, affecting the 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 [31]. Therefore, it would be advisable to carry out more experiments of this type in the future to help refine these initial results.