Acoustic Source Term Interpolation in a Hybrid Aeroacoustic Simulation ofthe Human Phonation
Michael WEITZ1 , Stefan SCHODER1 , Paul MAURERLEHNER1 , Sebastian FALK2 , Michael DÖLLINGER2 , and Manfred
KALTENBACHER1(1) Institute of Mechanics and Mechatronics, TU Wien, Austria, michael.weitz@tuwien.ac.at(2) Department of Phoniatrics and Pediatric Audiology at the Department of Otorhinolaryngology Head & Neck Surgery,

University Hospital Erlangen, Germany, sebastian.falk@uk-erlangen.deAbstract
Voice research is carried out experimentally by applying physical replicas and numerical models. Experimental
investigations involve high personnel, material and thus financial costs. Furthermore, the scientific outcome islimited to the few parameters that are measured at specific positions in the larynx replica. In contrast to that,numerical models provide a much higher amount of results since parameters from all possible locations in themodel can be obtained and analyzed. We follow a hybrid aeroacoustic approach, which performs in a first step
an incompressible flow simulation on a computational grid which is capable of resolving all relevant turbulentscales. In a second step, we compute the acoustic source terms on the flow grid and perform a conservative interpolation to the acoustic grid on which we solve the perturbed convective wave equation to obtain the acoustic
field. Finally, we validate the simulation results with results obtained from measurements.
Keywords: Finite element method, Computational aeroacoustics, Human phonation

1INTRODUCTION

The aim of the present study is to validate the aeroacoustic workflow of the human phonation process regarding
a potential clinical application of an efficient computational model. This demand for computational efficiency
is taken into account by the hybrid aeroacoustic approach. The hybrid aeroacoustic approach separates the flowand acoustic computations. Thereby, the flow results act as aeroacoustic sources in the acoustic computations,which means that the flow couples only into the acoustic field. The aeroacoustic methods are implemented in
the software package CFS++ [4].

2SIMULATION MODEL

The hybrid aeroacoustic approach is a powerful method in aeroacoustics to decouple the flow and acousticcomputations. We apply this approach to human voice research in order to obtain a computationally efficientand robust sound prediction. In a first step an incompressible flow field, which resolves all relevant turbulent
scales, is simulated. The incompressible flow simulation and its validation of the CFD simulation is publishedin [12]. In a second step, we compute the acoustic source terms on the flow grid and perform a conservativeinterpolation [5] to the acoustic grid, on which we solve the perturbed convective wave equation (PCWE) toobtain the acoustic field.
Section 2.1 describes the investigated geometry and how the acoustic simulation domain is derived from theexperimental setup and the flow simulation geometry to comply with the hybrid workflow. The aeroacoustic
model and the numerical setup of the acoustic simulation is discussed in Sec. 2.2. Additionally, we focus onefficient numerical computation via non-conforming interfaces [6].

1780

2.1 Geometry
The simulated geometry is derived from the geometry of the experiments and follows the rules of the hybrid
aeroacoustic approach. Section 3.1 describes the experimental setup in detail [7]. Figure 1 describes the test
rig, with a rectangular subglottal channel, the vocal folds, the ventricular folds, and the rectangular supraglottal
channel. Microphones measured the sound radiation of the configuration at 1 m distance from the end of thesupraglottal channel.

CAA
propagation
regionventricularfolds

yz

flow

1000

subglottal
channelmicrophonesupraglottal channel
1000

x CFD region = CAA source region

19000
10vocal folds

microphone

microphone

Figure 1. Schematic of the larynx model: subglottal, glottal, supraglottal regions, and the sound propagation
region with microphone positions. The CFD (bluish) and the CAA (reddish) geometry is adapted to the experimental setup. All dimensions are in mm. Extracted and modified from [7].
The CFD domain is depicted by the bluish region in Fig. 1. It includes the rectangular subglottal channel,
the vocal folds, the ventricular folds, and the rectangular supraglottal channel, that is, all parts of the test rig
where the flow speed is sufficiently high. The incompressible aerodynamics of the present validation case wassimulated with Star-CCM+ (Siemens PLM software, Plano, TX/USA) and the movement of the vocal folds
is prescribed with a driving frequency of f = 148 Hz symmetrically about the center line between the vocalfolds. Additionally to the CFD domain, the CAA geometry (reddish region) models the sound radiation to themicrophone positions.
2.2 Aeroacoustic model
Hardin and Pope [2] introduced the acoustic/viscous splitting technique for the prediction of flow-induced sound.
Afterwards, many scientists applied this idea and derived linear and non-linear wave equations [14, 1, 13, 11].
The essence of all these methods is that the flow field quantities are split into compressible and incompressible
parts

1781

pv
ρp̄ + pic + pc = p̄ + pic + pa
= v̄v + v ic + v c = v̄v + v ic + v a(1)a(3)== ρ̄ + ρ1 + ρ .

(2)

2.2.1 Governing equation
In this sense, the field variables are decomposed into mean ( p̄, v̄v, ρ̄) and fluctuating parts. Additionally, thefluctuating parts are further split into acoustic (pa , v a , ρ a ) and flow components (pic , v ic ). Finally, a density
correction ρ1 is built according to Eq. (3). Introducing an arbitrary Lagrangian-Eulerian description for theoperators∂
D=+ v − v rg · ∇ ,(4)Dt∂twhere v rg is the relative velocity of the grid, we arrive at the PCWE for moving meshes (see [5])1 Dpic1 D2 ψ aa−∆ψ
=−.
c2 Dt 2ρ̄c2 Dt(5)

This scalar convective wave equation is computationally efficient and describes aeroacoustic sources generated
by incompressible flow structures and wave propagation through moving media. This formulation reduces theaaa
number of unknowns (acoustic pressure pa = ρ̄ Dψ
Dt and particle velocity v = −∇ψ ) to just one scalar unknown,athe acoustic velocity potential ψ .
2.3 Numerical setup
The results of the aerodynamics simulation, which are required for the source term computation, are exported
as Ensight Gold-files. Every tenth CFD time step is written to the files leading to a temporal resolution of thesource term of 10 µs. The incompressible fluid density is set to ρ = 1.204 kg/m3 and the speed of sound isgiven by c = 343 m/s.
2.3.1 Domain
Figure 2 shows the computational domain of the acoustic simulation. As a consequence of the hybrid workflow
and the non-conforming meshing strategy, the acoustic domain consists of three separate parts: the larynx, thesimple rectangular vocal tract, and the propagation region. The larynx and the vocal tract together match theaerodynamic domain and contain aeroacoustic source. Furthermore, we have to note again that we use the simple rectangular vocal tract for the validation purpose with dimensions taken from the experimental study [7].
The propagation domain, a rectangular box located at the end of the vocal tract, is necessary to calculate thesound radiation to the microphone position in a 45◦ angle and a distance of 80 mm to the center line. The
distance of the microphone in the simulation is intentionally chosen smaller than the distance in the measurements to make the computation even more efficient. To account for the discrepancy of distances, we scale theresulting acoustic pressure of the simulation with 1/r, where r is the distance from the end of the supraglottal
channel to the respective microphone position.
2.3.2 Discretization
The existing capabilities of the simulation software allow well suited meshes for the used subdomains. Inside
the larynx, a mesh consisting of tetrahedral elements approximates the continuous geometry by 266000 elementswith a maximum element size of 1 mm. The discretization is chosen considerably different to the CFD grid, but
the utilized conservative interpolation strategy allows a CFD grid independent meshing. Inside the glottis region,
the characteristic length of the acoustic mesh is in the range of the aerodynamic discretization. The minimum1782

time step of 10 µs is enforced by the data exchange between the flow and the acoustic simulation. Based on theestimation formula for characteristic acoustic element length haha =λc== 3.4 mm ,20 20 fmax

(6)

we resolve a maximum frequency of fmax = 5 kHz with 20 linear finite elements per acoustic wave length λ tominimize the dispersion error if the actual element size is below the characteristic acoustic element length.

PML

Propagation
domainGlottisVocaltract

Non-conforming interfaces

PML

Propagation domain

Sound hard boundary

Figure 2. Geometry and mesh of the acoustic simulation: subglottal, glottal, supraglottal region, the soundpropagation region, and the perfectly matched layer (PML) regions.
The simple vocal tract, the propagation region, and the perfectly matched layer (PML) regions are meshedwith hexahedron elements. Table 1 summarizes the number of elements and the maximum element size in thedifferent regions. The total number of finite elements of the acoustic mesh sums up to about 1 million.
For simpler mesh generation, the different meshes are mathematically connected by non-conforming interfacesof Nitsche-type mortaring. The individual mesh generation reduces the number of total elements significantly
while maintaining accurate meshes for the wave equation. Therefore, this meshing strategy allows us to lower
the computational effort in the acoustic simulation.
We neglect the movement of the vocal folds and the convective term of the aeroacoustic source since the averaged flow speed is relatively low [12]. These assumptions reduce the grid complexity and the data transferred
as well as the computational cost.

1783

Table 1. Mesh and element statistics of the subdomains.
RegionNumber of Elements

max. Element Size in mmGlottisSimple vocal tract
Propagation region
Propagation region subglottal
PML subglottal
PML propagation
Total

2660001250
511700600600217750997900

1
3.6
3.6
1.5
1.5
3.6
3.6

2.3.3 Boundaries
All channel walls are modeled as fully reflecting, which is justified by a sufficiently high acoustic impedance
jump between air and tissue. Two PML regions ensure free field radiation without reflection [3]; one PML is
located at the inflow in the subglottal channel and the other PML surrounds the propagation region (see Fig.
2). The PML and the propagation region are connected by a conforming interface. It should be noted thoughthat the lower side of the propagation region is modeled as acoustically hard wall without PML.

3VALIDATION

3.1 Validation data
The experimental setup for generating the validation sound signal includes a mass flow generator with a silencer
installed behind it, the subglottal channel representing the synthetic trachea, the mounting device for the synthetic vocal folds model, and the supraglottal channel including false vocal folds [9]. The experimental setup
is designed in human length-scale and produces aerodynamically driven vocal fold vibrations as it is the caseduring normal speech.
The oscillations of the synthetic vocal fold models are flow-induced, and the vocal folds oscillate symmetrically with a fundamental frequency of f0 = 148 Hz. The entire larynx model has a rectangular cross-sectionof 15 mm × 18 mm, and its outlet of the supraglottal channel is located in a distance of 190 mm from the vocal folds. The experimental synthetic larynx model setup was multiply used and validated in previous studies[9, 7, 10, 8].
The experimental measurements of the sound signal were carried out in an anechoic room of the Chair ofSensor Technologies of the University Erlangen-Nürnberg. Four Brüel & Kjær 4189-L-001 (Naerum, Denmark)1/2 ” microphones with spherical characters were used to measure the sound pressure. The microphones werepositioned in the far field arranged around the channel exit and inclined by 45◦ concerning the channel axis asdisplayed in Fig. 3. The distance to the outlet of the supraglottal channel was 1 m (see Fig. 3b). The acousticsignal was amplified by a Brüel & Kjær Nexus (Naerum, Denmark) conditioning amplifier and finally A/Dconverted by a National Instruments AC-NI PXIe-4496 (Austin, Texas) acquisition board with a resolution of24 bits. All acoustic signals were sampled for 60 s using a sampling frequency fs = 96 kHz.
For the validation of the acoustic simulation, the experimental sound data were taken from [9]. More details
about the experimental setup and measurements can be found in [9, 7].
3.2 Results and Discussion
A comparison of the amplitude spectral densities (ASD) of both the simulations and measurements is depicted
in Fig. 4. The base frequency f0 and its harmonics are very well represented in the simulations. In the lower
frequency range, not only the respective frequencies are in good agreement with the measurement results but1784

(a) Top view.

(b) Side view.

Figure 3. Schematic of the experimental setup for sound measurements of vibrating vocal folds. All dimensions
in m. Extracted from [7].
also the amplitudes. In the higher frequency range, the simulation results start to show deviations from themeasurements, but the frequency of the majority of harmonics of the base frequency f0 is still acceptable.
Overall, the presented results show the efficiency and accuracy of the hybrid aeroacoustic workflow applied tothe acoustic simulation of the human phonation.

4CONCLUSIONSIn this paper, we have shown how a hybrid aeroacoustic workflow can be established for the human phonation process to ensure a highly efficient simulation and accurate results. Based on CFD results obtained on a
restricted subdomain, the source terms have been calculated for the perturbed convective wave equation. In anext step, these source terms have been interpolated to a mesh specifically suited for the acoustic propagation
simulation using a sophisticated conservative source term interpolation technique. Eventually, the results of theacoustic propagation simulation have been compared with validation measurements. Thereby, it was shown thatthe proposed simulation chain leads to accurate results.

1785

Figure 4. Comparison of the amplitude spectral density of the simulation results (blue) and measurement results(orange).

1786

ACKNOWLEDGEMENTS
The authors acknowledge support from the German Research Foundation (DFG) under DO 1247/10-1 (no.
391215328) and the Austrian Research Council (FWF) under no. I 3702.

REFERENCES[1] R. Ewert and W. Schröder. Acoustic perturbation equations based on flow decomposition via source filtering. Journal of Computational Physics, 188(2):365–398, 2003.
[2] J. Hardin and D. Pope. An acoustic/viscous splitting technique for computational aeroacoustics. Theoreticaland Computational Fluid Dynamics, 6(5-6):323–340, 1994.
[3] B. Kaltenbacher, M. Kaltenbacher, and I. Sim. A modified and stable version of a perfectly matched layertechnique for the 3-d second order wave equation in time domain with an application to aeroacoustics.
Journal of computational physics, 235:407–422, 2013.
[4] M. Kaltenbacher. Numerical Simulation of Mechatronic Sensors and Actuators: Finite Elements for Computational Multiphysics. Springer Berlin Heidelberg, 2015.
[5] M. Kaltenbacher. Computational Acoustics. CISM International Centre for Mechanical Sciences. Springer
International Publishing, 2017.
[6] M. Kaltenbacher, A. Hüppe, J. Grabinger, and B. Wohlmuth. Modeling and Finite Element Formulation
for Acoustic Problems Including Rotating Domains. AIAA Journal, 2016.
[7] S. Kniesburges. Fluid-Structure-Acoustic Interaction During Phonation in a Synthetic Larynx Model. Shaker,
2014.
[8] S. Kniesburges, V. Birk, A. Lodermeyer, A. Schützenberger, C. Bohr, and S. Becker. Effect of the ventricular folds in a synthetic larynx model. Journal of Biomechanics, 55:128–133, 2017.
[9] S. Kniesburges, A. Lodermeyer, S. Becker, M. Traxdorf, and M. Döllinger. The mechanisms of subharmonic tone generation in a synthetic larynx model. The Journal of the Acoustical Society of America,
139(6):3182–3192, 2016.
[10] A. Lodermeyer, S. Becker, M. Döllinger, and S. Kniesburges. Phase-locked flow field analysis in a synthetic human larynx model. Experiments in Fluids, 56(4):77, 2015.
[11] M. Munz, C.-D.and Dumbser and S. Roller. Linearized acoustic perturbation equations for low Mach
number flow with variable density and temperature. Journal of Computational Physics, 224(1):352–364,2007.
[12] H. Sadeghi, S. Kniesburges, M. Kaltenbacher, A. Schützenberger, and M. Döllinger. Computational modelsof laryngeal aerodynamics: Potentials and numerical costs. Journal of Voice, 2018.
[13] J. Seo and Y. J. Moon. Perturbed compressible equations for aeroacoustic noise prediction at low machnumbers. AIAA journal, 43(8):1716–1724, 2005.
[14] W. Z. Shen and J. N. Sørensen. Aeroacoustic modelling of low-speed flows. Theoretical and Computational
Fluid Dynamics, 13(4):271–289, 1999.

1787