Strojniški vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372 © 2014 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2014.1837	Special Issue, Original Scientific Paper
Received for review: 2013-12-13 Received revised form: 2014-02-14 Accepted for publication: 2014-04-01
Finite Element Formulations Applied to Outer Ear Modeling
Gaia Volandri* - Costantino Carmignani - Francesca Di Puccio - Paola Forte
University of Pisa, Department of Civil and Industrial Engineering, Italy
The work described in this paper is part of a broader research activity on the development of a virtual ear. The present study focuses on the tympanic membrane and auditory canal modeling, which are important components in sound transmission. The standard finite element method (FEM) and an alternative method (the generalized FEM), suitable for modeling sound propagation at high frequencies, were applied. Two domains (fluid and structural) for the auditory canal and the tympanic membrane, respectively, were considered in order to evaluate the coupling of the different methods and to apply a fluid-structure interaction formulation. ANSYS® software was used for solving FEM analyses, while GFEM simulations were obtained by implementing the method in Wolfram Mathematica®. Simulation results include modal response, pressure distribution in the auditory canal and displacement distribution in the tympanic membrane. The identified modal frequencies of the auditory canal agree with published data reported in the literature. The validation of such method with standard FEM simulation at increasing mesh density shows that FEM is more suitable for simulations of the human ear in the audible frequency range, although the generalized formulation could be convenient if an ear model including the whole head or the ultrasound frequency range were investigated. Keywords: finite element method, auditory canal, simulation, sound transmission
0 INTRODUCTION
The work described in this paper is part of a broader research activity on the development of a model of the human hearing perception, a kind of "virtual ear". It deals with the analysis and simulation of the vibratory behavior of the auditory apparatus in the conventionally considered audible frequency range, 20 Hz to 20 kHz. In particular, the present study is focused on the auditory canal (AC) including the tympanic membrane (TM) that represents a fundamental portion of the "normal acoustic path" formed by the outer, middle and inner ear.
In the last decades, many bioengineering methods have been applied to simulate the dynamic behavior of some parts of the ear, both distributed or lumped parameter methods, such as those based on electromechanical analogy or multi-body dynamics. However, for the simulation at low frequencies (<10 kHz) of sound propagation in the AC and for the mechanical-acoustic transmission through the TM, the finite element method (FEM) is the most frequently used approach [1] to [7].
The first finite element (FE) models of the TM appeared in the '70s [1]; since then, the FEM has been widely employed to model the ear structures due to its remarkable capability for analyzing complex geometries and the mechanical properties of anisotropic and inhomogeneous materials. Models including also, at least, the middle ear were developed by Kelly et al. [2], by Koike et al. [3], by Gan et al. [4] to [6] and Lee an Chen [7]. In the last years, hybrid FE and multi-body models of the TM and middle ear were also proposed by the present authors [8] and [9].
For simulating wave propagation with a standard FE model, i.e. with polynomial element shape functions, the mesh should respect the "rule of thumb", commonly accepted for many wave problems, that there should be at least 10 nodes per wavelength X
[10].	For sound transmission in air, X ranges from 15.6 mm (20 kHz) to 15.6 m (20 Hz), approximately. This means that at high frequencies, conventional FEM modeling can have a high computational cost, as it requires a very dense discretization of the domain.
Except for some attempts to apply domain decomposition with parallel processing, alternative methods or generalizations of conventional FEM
[11]	have been proposed in the literature for low wavelength acoustic problems. These methods have high accuracy and a minor need of mesh refinement with respect to standard FEM. Among these methods there are the spectral methods, as the spectral element method (SEM) or spectral finite elements (SFE), based on the fast Fourier transform (FFT) [12] or on orthogonal polynomials [13], and the wave element methods (WEM), as the ultra-weak variational formulation method (UWVF) [14], the partition of unity method (PU(FE)M) [15] and the generalized finite element method (GFEM) [16] and [17].
In this paper a conventional FEM analysis, carried out on an approximated model of the system formed by the AC and the TM, is compared with the results of an alternative method: the GFEM was selected as suitable for further extension of modeling and simulation of three-dimensional sound propagation problems to higher frequencies or higher dimensions of problem domain. This method was implemented in a commercial code (Wolfram Mathematica®) and applied preliminarily to simplified geometries
*Corr. Author's Address: Department of Civil and Industrial Engineering, Largo Lazzarino 56122, Pisa, Italy, gaia.volandri@ing.unipi.it
363
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372
that approximate the anatomy of the outer ear for a comparison with standard FEM in ANSYS®.
1 OVERVIEW OF COMPARED FE METHODS
The present investigation is focused on different FE formulations compared in terms both of accuracy and computational time. In this section, the theory of the adopted standard and generalized formulations is reported in brief.
1.1 Theoretical Background of the Finite Element Methods
As it is well-known, the basic idea of the FEM is to divide the continuous domain of a problem in a discrete set of elementary subdomains (elements) where the field function (e.g. displacement) can be approximated by means of simple basis functions [11]. Such basis functions are typically continuous, with piecewise continuous derivatives, and moreover can be easily integrated. In order to guarantee the global continuity of the displacement field, adjacent elements should have same values along their boundaries, thus the basis functions also support an interpolatory solution.
Accordingly, the displacement u at a point x within an element e can be written as:
ue (x, t ) = N(x) a (t ),
(3)
ue (x) = £ Nk (x) c
k'
(1)
where Nk(x) are the basis functions and cek the interpolating coefficients.
Differences between the FE methods selected for comparison can be attributed mainly to the basis functions adopted in element formulations. The peculiar characteristics of standard FEM, though well known, are here reported to ease the comparison with the other method.
1.1.1 Standard Finite Element Method
In standard FEM, interpolating coefficients are the element nodal displacements a.ek (k is the node numbers), so that Eq. (1) becomes:
ue (x) = £ N (x) ak
(2)
and the shape functions are typically low (first or second) order polynomials. A more compact expression can be obtained introducing the matrix form of the basis functions N(x):
which introduces also the generalization to time-dependent problems.
It can be worth noting that in case of isoparametric elements the same basis functions are used for mapping the elements from a reference domain into the physical domain.
Accordingly, the element mass and stiffness matrices are obtained from the following integrals on the element domain Qe:
Me = J NTpNdD, K" = J BtDBdD, (4)
where p is the material density and B = SN, S being a differential operator for calculating strain (£e = Sue = SNae = Bae), and D an elasticity matrix [11]. FE codes calculate such integrals numerically, i.e.
J f(x)dQ = Xf(X,)i
(5)
in particular by means of Gaussian quadrature [11], which for polynomials gives the exact solution.
Various procedures exist for the refinement of finite element solutions: the local approximation can be improved by polynomials of increasingly higher degree (p version), or, having fixed the polynomial degree p (typically p < 2), by decreasing the mesh size h (h version) or increasing the mesh size h and the degree p of polynomials (hp-version) [11].
1.1.2 Generalized Finite Element Method
The generalized finite element method, GFEM, is a combination of the standard FEM and the partition of unity method (PU(FE)M), aimed at introducing additional terms in the approximating function which enhance the global behaviour of the solution, also reflecting the known information about the boundary value problem [16]. Thus the standard polynomial FE solution is enriched with special functions u(x) (or handbook functions) usually non-polynomials, that means:
u(x, t) = N(x)a(t) + u(x).	(6)
However such special functions must respect the global regularity constraints, and thus are obtained by some specific (solution or problem-dependent) enrichment or handbook function hf(x), multiplied by the partition of unity (PU) NPU(x), i.e.
e
k
k
364
Volandri, G. - Carmignani, C. - Di Puccio, F. - Forte, P.
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372
u(x, t ) = N (x)a(t ) + NPU (x)hf(x)a* (t )
(7)
where a* are nodal unknown parameters that adjust the enrichment. Since this concept is at the basis of the extended method, it is worth reminding that a partition of unity in a domain Q is a set of functions such that:
) = 1, Vx eQ.
(8)
Consequently any function fx) can be reproduced by its product with the functions ^(x).
Frequently, but not necessarily, the NPU functions correspond to the standard FE shape functions N, often linear or bilinear, with a "hat" shape, defined on patches (much wider than elements) and zero everywhere else. A 2D representation of patches and elements is shown in Fig. 1. The generalization to the 3D tetrahedral case is direct.
with dl = l— and l = 1,2,...,p, and similarly
P
, 2n = m-,
q
m = 1,2

and employed as
Fig. 1. Elements and patches
handbook functions hf(x) in Eq. (7).
Such a direct generalization of the 2D case provides a distribution of directions which are concentrated around the poles of an imaginary sphere. As stated in [18], it is impossible to obtain an equally spaced distribution of directions. However, the choice of the logic of distribution of directions represents an important issue since a different logic of selection of propagation directions, based on a priori knowledge about the solution, can allow reducing the number of directions required to obtain a given level of accuracy.
A main issue in the implementation of GFEM is the possible linear (or almost linear) relation between the added handbook functions and the standard basis FE ones with consequent ill-conditioning problems [16]. As reported by many authors, an increasing number of wave directions involves a higher result accuracy with the drawback of introducing ill-conditioned system matrices and requiring dedicated integration techniques. Thus, it is often necessary to introduce some checks in the implementation on the condition number of the matrices and, if needed, to update the direction number assigned to each node.
1.2 Fluid-Structure Interaction Formulation
The local approximability is increased due to the handbook functions, while maintaining the existing infrastructures of the FE codes. In fact, the essential boundary conditions can be imposed in GFEM exactly as in standard FEM.
In [17], GFEM is applied to the solution of the problem of Helmholtz with the evaluation of special plane wave and wave band functions and functions of Vekua as handbook functions and the conclusion is that the use of plane wave functions involves a lower computational burden without substantial variations in the asymptotic accuracy, compared to more complex functions.
In the 2D case the plane waves can be expressed as shown in [17]. For 3D problems, as in the simulation of the fluid domain corresponding to the auditory canal, an extension is required and the plane waves can be expressed as a direct generalization of the 2D case, detailed in [18]:
s(i) _ gik(xsinO/ cos^m + y smQj sin$m + zcosOj )
(9)
The coupling between partial domains of the whole model represents a significant aspect of modeling that includes the debated issue of the fluid structure interaction (FSI) formulation. The fluid structure interaction at the domain interface implies that the acoustic pressure exerts a load on the structure and that the structural motion produces an effective load on the fluid. The introduction, in the system dynamics governing equations, of a coupling matrix accounting for the effective surface area and the normal to the interface area, represents a possible FSI formulation [19].
In the present study, the FSI coupling formulation described in [20] was adopted at the interface between fluid and structural domains. Such a formulation involves the building of an asymmetric element matrix for the interface elements (typically fluid elements), which have the pressure and three translational DoFs. For the interface fluid element, the structural mass, damping and stiffness matrices, as well as the fluid mass, damping and stiffness matrices,
Finite Element Formulations Applied to Outer Ear Modeling
365
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372
assume the typical FE form. Mass (MFS) and stiffness (Kfs) coupling matrices are defined according to the following equations:
Mfs =p0 At , Kfs = -A,	(10)
where p0 is the fluid density, and the coupling matrix A is obtained by the following integration on the interface surface S:
A = J NnT N'TdS,	(11)
S
where N and N' represent the fluid and structural shape function matrices, respectively, and the normal to the interface is indicated with n.
2 OUTER EAR MODELS
2.1 Anatomy
The auditory canal belongs to the outer portion of the ear, the tympanic membrane is instead generally included in what is called the middle ear, being at the interface of the outer and middle ear (Fig. 2).
Fig. 2. Schematic drawing of the human ear
The AC conveys the vibratory waves propagating in air to the middle ear. The morphology of the AC, though often approximated with a cylindrical geometry, typically presents two curves [21]; this shape (referred to as "S") not only facilitates the channeling of the wave but introduces variations (typically amplifications) at some resonance frequencies. Although inter-subject biological differences exist, there is a general agreement on the value of adult AC length of 25 to 32 mm. The cross-sectional area ranges from 65.45 to 75.53 mm2 at the TM to 90.13 to 96.16 mm2 at the canal entrance [4] and [7].
The understanding of the peculiarities of the AC that most influence the transmission of the signal and its coupling to the middle ear are important in the design of prostheses and in reconstructive surgery.
The TM, located at the end side of the ear canal, forms an angle of about 140° with the upper and lower walls of the channel; such an orientation gives a useful area of about 85 mm2, greater than the orthogonal cross-section of the AC.
The TM has a typical conical shape with an opening angle of 132 to 137° and a height of the cone of about 1.42 to 2 mm. The elliptical base of the cone has a vertical axis length ranging from 8.5 to 10 mm and a horizontal axis length ranging from 8 to 9 mm, with the apex (named umbo and assumed as reference point) facing the medial side [22].
The thickness of the TM is a critical parameter for modeling, since to date accurate detailed experimental measurements of the thickness distribution are not available and since it presents a high inter-subject variability in terms of absolute values. Although there is a general agreement in estimating a nonuniform thickness of the membrane, in modeling an approximated single thickness value, ranging from 30 to 150 ^m (with an average value of 74 ^m) for the human TM is usually adopted [4], [7] and [22].
The tissue of the TM is made of a multilayer structure with fibers oriented mainly in the radial and circumferential directions. Two main regions are usually distinguished a "Pars Tensa" (PT) and a "Pars Flaccida" (PF), having different size and mechanical properties.
The TM is connected on the medial side to the ossicular chain, while it is anchored along its circumference to the wall of the tympanic cavity by means of a fibro-cartilaginous ring (annular ligament).
2.2 Simplified Models of the Auditory Canal and Tympanic Membrane
2.2.1 Geometry and Material Properties of the AC
Two simplified models of the auditory canal were used that approximate the anatomy of an ear canal. A 22 mm long cylindrical and a pseudo-anatomical geometries (shown in Fig. 3) of a human auditory canal were adopted and imported into ANSYS® environment for the FEM analysis. The pseudo-anatomical geometry (including an air portion in the auricle) was extracted, through a semi-automatic segmentation algorithm, from computed tomography data provided by the Department of Otolaryngology II of Cisanello Hospital in Pisa. Size and morphology of the AC reconstructed geometry were in the literature range for adult healthy subjects (length of about 28 mm, average diameter of about 9 mm) [4].
366
Volandri, G. - Carmignani, C. - Di Puccio, F. - Forte, P.
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372
As regards the material properties, for the fluid contained in the AC the air medium was assumed as a compressible, inviscid fluid with uniform mean pressure and density; a 1.225 kg/m3 density value and a 340 m/s speed of sound value were, thus, adopted. The damping was not considered in this evaluation study.
The AC bony wall was simulated with clamped boundary conditions and a distributed uniform harmonic sound pressure load of 90 dB SPL (corresponding to 0.632 Pa) was applied at the inlet of the AC while the outlet was not constrained.
An isotropic material model with Young's modulus of 20 MPa was considered [1]. The Poisson's ratio was set equal to 0.3; density was assumed equal to 1.2x103 kg/m3.
Triangular shell elements with three nodes were used for the TM, when included.
As concerns the boundary conditions, the TM was fully clamped at the periphery, i.e. the annular ligament was not considered, as well as the ossicular chain.
2.3 Implemented Methods
2.2.2 Mesh Definition for the AC Model
Tetrahedral elements were employed for the cylindrical and the pseudo-anatomical models (Fig. 3), as it is typically adopted for complex biological geometries.
The size of the elements was set by means of a convergence criterion based on the value of the first natural frequencies up to 20 kHz; the percent relative error vs. the logarithm of number of DoFs was estimated (variation less than 1%) [11].
a)
b)
Fig. 3. a) Cylindrical and b) anatomical geometries and mesh of a human auditory canal
2.2.3 Tympanic Membrane Model
The TM was included in the model with the aim of accounting for the fluid-structure interaction.
Although the TM has a peculiar shape (See Section 2.1), a flat circular geometry was adopted with the aim of facilitating the definition of the fluid-structure interface. A 0.074 mm thickness was adopted, deduced from [4].
The implementation and comparison of the methods was carried out in Mathematica® environment on the AC approximated geometries.
Firstly, the GFEM implementation required the definition of elements and patches and the setting of the wavenumber.
The partition of unity and the shape functions Nk were chosen as the linear shape functions for the brick element with eight nodes of the standard FEM.
In this paper special plane wave functions, Wj, were implemented, whose more general expression is:
Wji) = e*",
(12)
where r is the position vector and k the wave vector.
The unit vectors of the propagation directions of the plane-wave functions, in the 3D problems, were selected based on a priori knowledge of the solution or following the optimized "spherical covering" logic, borrowed from another method (e.g. ultra weak variational formulation, UWVF) [14]. The optimized "spherical covering" logic identifies n points on an imaginary sphere, centered in each patch node, so that the maximum distance of each node belonging to the sphere from the nearest of the n considered points is minimized. Such a criterion allows obtaining a direction distribution as homogeneous as possible. In the present study, n was chosen in the 4 to 124 range due to requirements of computational cost and the same number of directions was assigned to each vertex of the patches, although the method requires neither an equally spaced distribution of directions nor an equal number of DoFs at each node.
The GFEM method was implemented for the simulation of the fluid domain while the FEM method was implemented, in Mathematica® as well, for the simulation of the structural domain (limited to the TM), when included.
The ANSYS® FLUID30 and SHELL63 (Discrete Kirchoff element, DKT) formulations were adopted
Finite Element Formulations Applied to Outer Ear Modeling
367
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372
and implemented for the fluid and structural elements, respectively. The ANSYS® FLUID30 element formulation was selected as suitable for fluid/structure interaction problems and sound wave propagation applications. The element has a brick shape with eight corner nodes and four DoFs per node: three displacements (only at nodes on the interface) and pressure. The element adopts linear shape functions without extra shape functions, and a standard Gaussian quadrature for bricks (2*2*2 points).
For the comparison and validation of the generalized method, the standard FEM solution at different levels of mesh refinement obtained with the commercial code ANSYS®, was used as a reference.
2.4 Comparison Tests
2.4.1	Modal and Harmonic Analyses on the AC Models
Numerical free-free or constrained modal and harmonic analyses were carried out on the models, intending as constrained and free-free conditions when sound pressure load is applied or not, respectively, at the AC inlet. Pressure distribution in the AC was evaluated as main result of the harmonic analysis.
In addition to a check on eigenfrequencies in a specific frequency range ([20 Hz to 20 kHz]), the modal shape correlation was quantified by the modal assurance criterion (MAC) [23]. MAC results are usually represented as a square matrix, correlating the reference modes (FEM, test, theory etc.) to verification modes. It assumes a nearly zero value in the presence of incompatible mode shapes, a unitary value in case of perfect correlation and intermediate values in case of partial correlation.
2.4.2	FSI of Auditory Canal and Tympanic Membrane
Harmonic analyses were also carried out on the cylindrical geometry with circular membrane accounting for the fluid-structure interaction. The cylindrical geometry was preferred since it allows a simple definition of interfaces for the fluid-structure formulation.
The material properties and the boundary conditions of the AC, except for the outlet, were set as previously described.
The harmonic analysis on the implemented models provides the distribution of pressure and normal displacement in the fluid and structural domain, respectively (coupled by FSI).
3 RESULTS
3.1 Comparison of Methods on the AC Model
In order to compare the two finite element formulations, free-free and constrained modal analyses were performed on the simplified geometries (Fig. 3) of the fluid domain. Eigenfrequencies in the 20 Hz to 20 kHz frequency range (typically the first three in the pseudo-anatomical case) were compared and mode shapes correlated, using the MAC matrix.
The modal frequencies obtained in the free-free and constrained modal analysis with FEM and GFEM in the cylindrical and pseudo-anatomical geometries are reported and compared in Table 1.
Table 1. First three modal frequencies; the results are shown for a), b) cylindrical geometries; and for c), d) pseudo-anatomical geometries;a), c) for the free-free modal analysis; b), d) for the constrained modal analysis
	f [Hz]	1°	2°	3°
a)	FEM	7724	15489	23257
	GFEM	6648	13880	24097
	error [%]	14	10	-3.6
b)	FEM	3864	11601	19367
	GFEM	3670	11269	19832
	error [%]	5	2.8	-2.4
c)	FEM	6616	12022	15258
	GFEM	5587	11960	13345
	error [%]	15	0.5	12.5
d)	FEM	3961	11103	17119
	GFEM	3330	11229	17677
	error [%]	16	-1.1	-3.2
It is worth noting that the modal frequencies obtained with the pseudo-anatomical geometry are in agreement as order of magnitude with the three resonance frequencies in the audible range (3176, 9528 and 15880 Hz) inherent to the external auditory canal, depending on its morphology and length, reported in [24]. Moreover, in the literature the first natural frequency inherent to the auditory canal is often identified in the 3 to 4 kHz range [2] and [24] confirming the estimated value. Finally, the results obtained by the two numerical approaches appear to differ majorly at low frequency where standard FEM is considered reliable and convenient while they are comparable at higher frequencies where GFEM should become more advantageous.
The comparison of the mode shapes was carried out in Mathematica® environment by importing the results obtained in the different software environments for the proposed methodologies. Fig. 4 shows, for the cylindrical (a, b) and for the pseudo-anatomical
368
Volandri, G. - Carmignani, C. - Di Puccio, F. - Forte, P.
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372
geometry (c, d), for the free-free (a, c) and constrained (b, d) modal analysis respectively, the matrix representations of the modal shape correlation of the two methods, for an equal number of nodes.
The maximum/minimum values on the MAC matrix diagonal as well as the out of diagonal maximum values are shown in Fig. 4 for the four above mentioned cases.
a)
fem 1 2 3
b)
0.90			0.85			
					0.08	
	0.18	0.60			0.66	
c) d)						
0 0.1 0.2 0.3 0,4 0,5 0,5 0.7 0.S 0.0 1
Fig. 4. Modal shape correlation by MAC; the results are shown for the a), b) cylindrical and c), d) pseudo-anatomical geometries, for
the free-free a), c) and constrained b), d) modal analysis
All cases show consistency in the first three natural frequencies obtained by the various methodologies, less in the mode shapes.
As an example of the harmonic analysis results, the pressure distribution at 10 kHz in the pseudo-anatomical model of the AC is shown in Fig. 5.
The GFEM results obtained with a coarse element mesh (144 nodes, 288 DoFs) are compared with the ANSYS® FEM results at equal and refined (3216 nodes, 3216 DoFs) mesh. The TM side is indicated in the figures. The GFEM and FEM harmonic results, in the audible frequency range agree within tight tolerances.
The GFEM is more expensive from the computational point of view, since GFEM presents the difficulty of calculation in the complex field and entails problems of ill-conditioning and linear dependence in the phase of construction of the element matrices, for which it often requires heavy integration techniques. However, with respect to alternative methods of the literature, GFEM presents
the advantage of keeping the standard FEM mesh, even with tetrahedral elements, with nodes belonging to the element boundary.



a)
b)
c)
Fig. 5. Harmonic analysis results at 10 kHz: pressure distribution in the pseudo-anatomical geometry with a) GFEM (coarse mesh) and b) coarse mesh, and c) FEM refined mesh
In conclusion, FEM is more suitable for outer ear simulations in the 20 Hz to 20 kHz frequency range. Moreover it is available in commercial codes.
However, further investigations elucidate that the GFEM can be more suitable for other applications involving higher frequencies (e.g. a 22 mm long cylindrical model at 100 kHz) or higher characteristic dimensions of the problem (e.g. a scaled/homothetic 22 cm long cylindrical model, at 10 kHz). The results of these further simulations are shown in Figs. 6 and 7, respectively, for different mesh density. The longitudinal section of the cylinder, instead of the cylindrical surface, is plotted in Figs. 6a and 7a.
As one can see, the results of Fig. 6a are comparable to those of Fig. 6c which represents the
Finite Element Formulations Applied to Outer Ear Modeling
369
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372
a)
c)
Fig. 6. Harmonic analysis in a cylindrical model (22 mm long) at 100 kHz with: a) GFEM (longitudinal section) and (b, c) FEM varying the mesh density: a) 137 nodes and 166 DoFs, b) 137 nodes and DoFs, c) 28487 nodes and DoFs
a)
c)
Fig. 7. Harmonic analysis in a cylindrical model (scaled 22 cm long cylinder) at 10 kHz with: a) GFEM (longitudinal section) and b), c) FEM varying the mesh density: a) 137 nodes and 166 DoFs, b) 137 nodes and DoFs, c) 28487 nodes and DoFs
FEM reference solution, obtained with a considerably higher number of degrees of freedom. Then, if Fig. 6a is compared with Fig. 6b, obtained with an equal number of degrees of freedom, the differences and the enhancement are evident. The same observations hold for Fig. 7.
The use of the implemented generalized element formulations can be convenient, especially as frequency increases (and therefore the wave number) because it allows to achieve, with a coarse mesh (which does not satisfy the rule of thumb of ten nodes per wavelength conventionally accepted for standard FEM), an accuracy comparable to that obtained with a fine mesh in standard FEM formulations.
These first indications suggest that the integration of these advanced techniques in a FE model of the ear could be useful if, for example, the whole head were included or if the ultra-sound field were investigated.
370
3.2 Comparison of Methods on the FSI Problem
Concerning the simulation of the fluid- structure problem by the combination of techniques GFEM and FEM for the fluid and structural domains, respectively, the distributions of pressure inside the AC and TM displacement at 200 Hz are compared with the results obtained with standard FEM in ANSYS® with a coarse and refined mesh, in Figs. 8 and 9, respectively.
As highlighted in the previous section, the results indicate that GFEM and FEM results, in the investigated frequency range, are in agreement within acceptable tolerances.
4 CONCLUSIONS
In this work some modeling aspects of the human ear canal and tympanic membrane were examined also
Volandri, G. - Carmignani, C. - Di Puccio, F. - Forte, P.
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372
a)
c)
Fig. 8. Harmonic analysis results at 200 Hz in the FSI problem: pressure distribution in the auditory canal with: a) GFEM (coarse mesh, longitudinal section) and b) FEM (coarse mesh), c) FEM (refined mesh
a)
c)
Fig. 9. Harmonic analysis results at 200 Hz in the FSI problem: normal displacement distribution in the tympanic membrane with: a) GFEM and b), FEM (coarse mesh), c) FEM (refined mesh)
considering the fluid-structural coupling that occurs between these two components of the acoustic path. Standard and generalized finite element models were implemented, tested and compared.
GFEM and FEM modal and harmonic results in the 20 Hz to 20 kHz range agree within tight tolerances in all tested cases. Thus, FEM appears more suitable for simulations in the audible frequency range, assuming the typical ear size of a human being, due to its relatively limited computational burden and the availability of commercial codes.
However, these preliminary results show also that the generalized finite element formulation can be convenient in short-wave acoustic problems with the aim of simulating the auditory apparatus including the whole head or investigating the ultra-sound field, i.e. when the wave length is shorter than the characteristic dimension of the problem.
The analysis results concerning the natural frequencies of the auditory canal are consistent with some published studies [2] and [24] which identify three modal frequencies in the audible range with the first modal one included in the 3 to 4 kHz range.
As future developments, the implementation and integration of this method in a commercial code can be planned. The application of GFEM to a complete ear model including bone conduction in the head can be considered. Moreover GFEM can be proposed for convenient application to other fields (e.g. ultrasounds in medicine).
5 ACKNOWLEDGMENTS
The support and contribution of Prof. Stefano Berrettini and Dr. Luca Bruschini of the U.O.
371
Finite Element Formulations Applied to Outer Ear Modeling
Strojniski vestnik - Journal of Mechanical Engineering 60(2014)5, 363-372
Otorinolaringoiatria 2° of the Cisanello hospital in
Pisa, Italy are gratefully acknowledged.
6 REFERENCES
[1]	Funnell, W.R.J., Laszlo, C. (1978). Modeling of the cat eardrum as a thin shell using the finite-element method. Journal of the Acoustical Society of America, vol. 63, no. 5, p. 1461-1467, D0I:10.1121/1.381892.
[2]	Kelly, D.J., Prendergast, P.J., Blayney, A.W. (2003). The effect of prosthesis design on vibration of the reconstructed ossicular chain: a comparative finite element analysis of four prostheses. Otology & Neurotology, vol. 24, no. 1, p. 11-19, D01:10.1097/00129492-200301000-00004.
[3]	Koike, T., Wada, H., Kobayashi, T. (2002), Modeling of the human middle ear using the finite-element method. Journal of the Acoustical Society ofAmerica, vol. 111, no. 3, p. 1306-1317, D01:10.1121/1.1451073.
[4]	Gan, R.Z., Feng, B., Sun, Q. (2004). Three-dimensional finite element modeling of human ear for sound transmission. Annals of Biomedical Engineering, vol. 32, no. 6, p. 847-859, D01:10.1023/ B:ABME.0000030260.22737.53.
[5]	Gan, R.Z., Sun, Q., Feng, B., Wood, M.W. (2006). Acoustic-structural coupled finite element analysis for sound transmission in human ear— Pressure distributions. Medical Engineering & Physics, vol. 28, no. 5, p. 395-404, D01:10.1016/j. medengphy.2005.07.018.
[6]	Gan, R.Z., Reeves, B.P., Wang, X. (2007). Modeling of sound transmission from ear canal to Cochlea. Annals of Biomedical Engineering, vol. 35, no. 12, p. 21802195, D0I:10.1007/s10439-007-9366-y.
[7]	Lee, C.F., Chen, P.R., Lee.W.J., Chou, Y.F., Chen, J.H., Liu, T.C. (2010). Computer aided modeling of human mastoid cavity biomechanics using finite element analysis. EURASIP Journal on Advances in Signal Processing, paper: 203037, D0I:10.1155/2010/203037.
[8]	Volandri, G., Di Puccio, F., Forte, P. (2012). A sensitivity study on a hybrid FE/MB human middle ear model. ASME 2012 11th Biennial Conference on Engineering Systems Design and Analysis, vol. 4, p. 217-226, D0I:10.1115/ESDA2012-82326.
[9]	Volandri, G., Di Puccio, F., Forte, P., Manetti, S. (2012). Model-oriented review and multi-body simulation of the ossicular chain of the human middle ear. Medical Engineering & Physics, vol. 34, no. 9, p. 1339-1355, D0I:10.1016/j.medengphy .2012.02.011.
[10]	Zienkiewicz, O.C., Taylor, R.L., Nithiarasu, P. (2005). The Finite Element Method for Fluid Dynamics. Elsvier Butterworth-Heinemann, Burlington.
[11]	Zienkiewicz, O.C., Taylor, R.L. (2000). The Finite Element Method (5th ed.) Volume 1 - The basis, Elsevier Butterworth Heinemann, Burlington.
[12]	Doyle, J.F. (1997). Wave propagation in Structures, 2nd ed., Springer Science+Business Media New York .
[13]	Peng, H., Meng, G., Li, F. (2009). Modeling of wave propagation in plate structures using three- dimensional spectral element method for damage detection. Journal of Sound and Vibration, vol. 320, no. 4-6, p. 942-954, D0I:10.1016/j.jsv.2008.09.005.
[14]	Huttunen, T., Gamallo, P., Astley, R.J. (2009). Comparison of two wave element methods for the Helmholtz problem. Communications in Numerical Methods in Engineering, vol. 25, no. 1, p. 35-52, D01:10.1002/cnm.1102.
[15]	Gamallo, P., Astley, R.J. (2006). The partition of unity finite element method for short wave acoustic propagation on non-uniform potential flows. International Journal for Numerical Methods in Engineering, vol. 65, no. 3, p. 425-444, D01:10.1002/ nme.1459.
[16]	Strouboulis, T., Babuska, I., Copps, K. (2000). The design and analysis of the generalized finite element method. Computer Methods in Applied Mechanics and Engineering, vol. 181, no. 1-3, p. 43-69, D01:10.1016/ S0045-7825(99)00072-9.
[17]	Strouboulis, T., Hidajat, R., Babuska, I. (2008). The generalized finite element method for Helmholtz equation. Part II: Effect of choice of handbook functions, error due to absorbing boundary conditions and its assessment. Computer Methods in Applied Mechanics and Engineering, vol. 197, no. 5, p. 364380, D01:10.1016/j.cma.2007.05.019.
[18]	Laghrouche, O., Bettess, P., Perrey-Debain, E., Trevelyan, J. (2003). Plane wave basis finite-elements for wave scattering in three dimensions. Communications in Numerical Methods in Engineering, vol. 19, no. 9, p. 715-723, D01:10.1002/cnm.632.
[19]	Fluids Analysis Guide (2005). ANSYS Release 10.0, Ansys Inc., Canonsburg.
[20]	Theory Reference (2004). ANSYS Release 9.0, Ansys Inc., Canonsburg.
[21]	Gray, H. (2000). Anatomy of the Human Body. 20th edition, Lewis, W.H. (ed.) Lea and Febiger, Philadelphia & Bartleby.com, New York, from http:// www.bartleby.com/107/ accessed on 2013-01-01.
[22]	Volandri, G., Di Puccio, F., Forte, P., Carmignani, C. (2011). Biomechanics of the tympanic membrane. Journal of Biomechanics, vol. 44, no. 7, p. 1219-1236, D0I:10.1016/j.jbiomech.2010.12.023.
[23]	Ewins, D.J. (2000). Model validation: Correlation for updating. Sadhana, vol. 25, no. 3, p. 221-234, D0I:10.1016/j.jbiomech.2010.12.023.
[24]	Vallejo, L.A., Delgado, V.M., Hidalgo, A., Gil-Carcedo, E., Gil-Carcedo, L.M., Montoya, F. (2006). Modeling of the geometry of the external auditory canal by the finite elements method. Acta Otorrinolaringologica Espa-ola, vol. 57, no. 2, p. 82-89, D0I:10.1016/S0001-6519(06)78667-8. (in Spanish)
372
Volandri, G. - Carmignani, C. - Di Puccio, F. - Forte, P.