This article deals with large-eddy simulations of three-dimensional incompressible laryngeal flow followed by acoustic simulations of human phonation of five cardinal English vowels, /ɑ, æ, i, o, u/. The flow and aeroacoustic simulations were performed in OpenFOAM and in-house code openCFS, respectively. Given the large variety of scales in the flow and acoustics, the simulation is separated into two steps: (1) computing the flow in the larynx using the finite volume method on a fine moving grid with 2.2 million elements, followed by (2) computing the sound sources separately and wave propagation to the radiation zone around the mouth using the finite element method on a coarse static grid with 33 000 elements. The numerical results showed that the anisotropic minimum dissipation model, which is not well known since it is not available in common CFD software, predicted stronger sound pressure levels at higher harmonics, and especially at first two formants, than the wall-adapting local eddy-viscosity model. The model on turbulent flow in the larynx was employed and a positive impact on the quality of simulated vowels was found.

Aeroacoustic simulations undoubtedly have strong potential to be applied in clinical diagnostics, treatment control, and to support medical education. Regarding this application, computer analysis can provide highly resolved three-dimensional (3D) data of the flow and acoustic field for further studies. Such highly resolved 3D data are infeasible by in vivo, excised larynx or synthetic vocal fold measurements (Döllinger , 2011; Kniesburges , 2011).

Although numerical tools are available, realistic 3D turbulent flow simulations are computationally expensive as the time-consuming accurate modeling of supraglottal turbulence is essential for flow-induced sound generation (Mattheus and Brücker, 2011). Conventional turbulence modeling approaches, such as unsteady Reynolds-averaged Navier-Stokes (uRANS) equations are inadequate for aeroacoustics because they do not provide the instantaneous state of flow quantities. Therefore, in comparison with a direct numerical simulation (DNS), large-eddy simulation (LES) is a beneficial balance between the resolution of turbulent structures and the accurate prediction of turbulent sound generation mechanism. One of the first studies employing LES to research laryngeal aeroacoustics was the work of Suh and Frankel (2007), who combined a compressible LES and acoustic analogy published by Ffowcs Williams and Hawkings (1969; FW-H) in a static model of the human glottis for human voice signal predictions. Mihaescu (2010) employed the LES capability to study the laryngeal airflow during phonation and inspiration. Schwarze (2011) explored an implicit LES, where the intrinsic dissipation of the numerical method was assumed to act as a subgrid-scale model. During recent years, simulations have advanced considerably (Bodaghi , 2021; Tokuda and Shimamura, 2017; Yoshinaga , 2017) toward a state where full-scale aeroacoustic simulations on realistic computed tomography (CT)- or magnetic resonance imaging (MRI)-based geometries are possible. It can be anticipated that these simulations could be used for subject-specific presurgical predictions of vocal fold oscillations (Avhad , 2022). These phonation simulations can be useful and help to improve the voice quality for subjects suffering from various vocal fold dysfunctions (Falk , 2021; Sadeghi , 2019a,b) or evaluate potential effects on voice production affected by an implant insertion in medialization laryngoplasty (Zhang , 2020).

In clinical applications, the reliability of the LES is a key factor. Hence, further studies on the accuracy of subgrid-scale models are important. In our previous study (Lasota , 2021), the implementation of a relatively new subgrid-scale anisotropic minimum dissipation (AMD) model, published first by Rozema (2015), was described. The reasons for introducing the given model are based on the major features of the minimum dissipation models, which are constructed to confine time derivative of subgrid-scale energy and satisfy the need to dissipate the energy of subgrid scales. In addition, the AMD model is also derived for anisotropic meshes, which can circumvent an overuse of corrections needed to retain the stability of the numerical schemes when severe mesh element deformation occurs due to vocal fold motion. Rozema (2015) confirmed high consistency between results obtained by AMD and the DNS of turbulent channel flow and between AMD and the DNS of temporal mixing layers, where both were performed on anisotropic grids. These flow regimes are comparable to the flow structures arising inside the larynx.

In this work, AMD is investigated for biological laryngeal flows based on the above-mentioned advantages. The flow field prediction capabilities are assessed and compared to state-of-the-art LES using a conventional one-equation (OE) and wall-adapting local eddy-viscosity (WALE) subgrid-scale turbulence models. Furthermore, the AMD model's impact on the aeroacoustic sound generation and the produced voice signal is evaluated for the first time.

The article is organized into two major sections describing the computational fluid dynamic (CFD) and computational aeroacoustic (CAA) parts. CFD in Sec. II describes models, geometry, boundary conditions, mesh, discretization, numerical solution, and results. CAA in Sec. III keeps the same structure.

LES is a mathematical concept for modeling turbulent flows, which deals with flow structures carrying most kinetic energy, k, i.e., spatially organized large scales. These consist of two main categories: coherent structures and coherent vortices of recognizable shape (Lesieur , 2005). In the numerical implementation, the characteristic length, Δ, defining a cutoff between resolved large scales and modeled subgrid scales, is usually given by the mesh grid spacing (Versteeg and Malalasekera, 2007).

In the LES concept, any flow variable f ( x , t ), where x = ( x 1 , x 2 , x 3 ) is the spatial coordinate and t is time, may be decomposed as
f ( x , t ) = f ¯ ( x , t ) + f ( x , t ) ,
(1)
where f ¯ ( x , t ) = G f ( x ) f ( x , t ) = G f ( r , x , Δ ) f ( x r , t ) d r is the large-scale component obtained by spatial filtering, and f ( x , t ) is the small subgrid-scale contribution. Filtered variables for LES are functions of time and space, unlike the Reynolds-averaged variables, hence, in LES, f ¯ ¯ f ¯ , f ¯ 0. The convolution introduced above contains a filter function, Gf, separating spatial scales. The filter used in this study is the top-hat filter, which is a common choice in low-order finite volume methods,
G f ( r , x , Δ ) = { 1 / Δ 3 for | r | Δ / 2 , 0 , otherwise .
(2)
The continuity and momentum equations for the incompressible fluid flow with LES filtering applied can be written as
u ¯ i x i = 0 , u ¯ i t + x j ( u i u j ¯ ) = 1 ρ p ¯ x i + ν 2 u ¯ i x j x j ,
(3)
where u ¯ i is the filtered velocity, p ¯ represents the filtered static pressure, and ν is the kinematic molecular viscosity. The term u i u j ¯ is the dyadic product and cannot be expressed directly (Ferziger, 1998). Modification of the momentum equation (3) by + ( / x j ) ( u ¯ i u ¯ j ) yields
u ¯ i t + x j ( u ¯ i u ¯ j ) = 1 ρ p ¯ x i + ν 2 u ¯ i x j x j τ i j x j .
(4)
The new term on the right-hand-side of Eq. (4) is the divergence of the subgrid-scale turbulent stress tensor,
τ i j = u i u j ¯ u ¯ i u ¯ j ,
(5)
and left to be modeled to close the set of equations. Since the turbulence is not fully understood, a wide range of closure subgrid-scale models have been introduced, often using heuristic and ad hoc techniques.

1. OE subgrid-scale model

The OE model derived by Yoshizawa and Horiuti (1985) computes the transport equation for the turbulent kinetic subgrid-scale energy, kSGS, such that
k SGS t + u ¯ j k SGS x j x j ( ( ν + ν t ) k SGS x j ) = τ i j S ¯ i j C ϵ k SGS 3 / 2 Δ ,
(6)
where C ϵ = 1.048 is the model constant, S ¯ i j is the resolved rate-of-strain tensor (symmetric part of the velocity gradient tensor), and νt is the turbulent viscosity.
Unlike the Smagorinsky model (Smagorinsky, 1963), which disregards the first three terms in Eq. (6), OE also considers the history effects for kSGS. The production term, τ i j S ¯ i j, modeling the decay of turbulence from the resolved scales to the subgrid scales via the energy cascade, is approximated by
τ i j S ¯ i j = 2 ν t S ¯ i j S ¯ i j .
(7)
OE relies on the subgrid-scale eddy viscosity,
ν t O = C ν Δ k SGS ,
(8)
where C ν = 0.094 due to the Kolmogorov law.

Neither the Smagorinsky nor the OE model can reproduce the laminar to turbulent transition, and both tend to overpredict the production rate and, thus, the turbulent viscosity in free shear layers and near the walls. This is caused by the fact that the term S ¯ i j S ¯ i j is large in the regions of pure shear because it is only related to the rate-of-strain, S ¯ i j, and not to the rate-of-rotation, Ω ¯ i j (Lesieur , 2005).

2. WALE subgrid-scale model

The inaccuracy concerning free shear and boundary layer treatment, caused by the OE model, can be alleviated by using the WALE model (Nicoud and Ducros, 1999). WALE is able to model accurately turbulent structures with strain or rotation rate or even both. Thanks to this behavior, the pure shear flow located near solid boundaries will cause the eddy-viscosity to vanish (Nicoud and Ducros, 1999). The turbulent viscosity computed by WALE is defined as
ν t W = C k Δ k SGS W ,
(9)
where the turbulent energy of subgrid scales, k SGS W, is
k SGS W = ( C w 2 Δ C k ) 2 ( s i j d s i j d ) 3 ( ( S ¯ i j S ¯ i j ) 5 / 2 + ( s i j d s i j d ) 5 / 4 ) 2 .
(10)
The model constants are set to C w = 0.325 and C k = 0.094, and s i j d is the traceless symmetric part of the square of the velocity gradient. WALE accounts for the rate of rotation in the computation of Eq. (9) and, thus, the turbulent viscosity tends to zero near walls. Hence, it is not necessary to use any ad hoc damping methods.

3. AMD subgrid-scale model

The main objective of the model is to ensure that the energy of subgrid scales, kSGS, is not increasing,
t Ω Δ 1 2 u i u i d x 0.
(11)
AMD was derived in Rozema (2015) with modified Poincaré inequality addressing the grid anisotropy. If subgrid scales are assumed to be periodical on the filter box, Ω Δ, it is possible to apply the Poincaré inequality and, thus, define the upper bound of kSGS such that
Ω Δ 1 2 u i u i d x C Ω Δ 1 2 ( i u j ) ( i u j ) R 1 d x .
(12)
The term R1 in Eq. (12) corresponds to the velocity gradient energy, and C is the Poincaré constant, C = ( Δ / π ) 2, for the LES filter of width Δ.
The AMD model can sidestep the dependence of the model constant on Δ by using the modified Poincaré inequality,
Ω Δ 1 2 u i u i d x C A Ω Δ 1 2 ( Δ x i i R 3 u j ) ( Δ x i i u j ) R 2 d x ,
(13)
where the filter box, Ω Δ, has dimensions Δ x 1 , Δ x 2, and Δ x 3, and CA is the modified Poincaré model constant. The term R2 is the scaled velocity gradient energy and R3 is the scaled gradient operator. The inequality (13) demonstrates that the subgrid energy is confined by imposing a bound on the term R2. Time derivative is applied on the term R2, and the evolution equation of R2 on the filter box is expressed as
t ( 1 2 ( Δ x i i u j ) ( Δ x i i u j ) ) = ( Δ x k k u i ) ( Δ x k k u j ) S i j R 4 ( ν + ν t A ) Δ x k k ( i u j ) Δ x k k ( i u j ) + i f i ,
(14)
where the term R4 is the production of the scaled velocity gradient energy. The following inequality in Eq. (15) ensures that the AMD model predicts sufficient dissipation to stop the production of scaled velocity gradient energy, R4:
Ω Δ R 4 d x ν t A C A Ω Δ ( i u j ) ( i u j ) d x ,
(15)
where the minimum dissipation effect is ensured by satisfying
ν t A = C A max { Ω Δ ( Δ x k k u i ) ( Δ x k k u j ) S i j d x , 0 } Ω Δ ( l u m ) ( l u m ) d x .
(16)
Integrals in Eq. (16) can be approximated by the mid-point rule, and the turbulent viscosity for AMD, ν t A, results in a more practical form,
ν t A = C A max { ( Δ x k k u i ) ( Δ x k k u j ) S i j R 4 , 0 } ( l u m ) ( l u m ) R 5 ,
(17)
where the terms in vector notation are
R 4 = ( Δ x u ) ( Δ x u ) : S , R 5 = ( u ) : ( u ) .
(18)

The AMD model was implemented into the OpenFOAM code as the LAAMD library.1

Based on decaying grid turbulence tests, the value of the constant, CA, in Eq. (17) is recommended with respect to the order of discretization of Navier-Stokes equations, which were tested on decaying grid turbulence cases. Rozema (2015) has concluded that AMD gave the best results with CA = 0.3 for a central second-order scheme and C A = 0.212 for a fourth-order scheme. A recent study by Zahiri and Roohi (2019) states an optimal value of the constant C A = 1 / 3 = 0.577 based on various test cases. Lasota (2022) performed his own tests with model constants, C A = 0.3 and C A = 0.57735, on cases with turbulent flow in a planar channel and periodic hill. Based on better agreement with velocity profiles obtained by DNS or experiments, C A = 0.3 was chosen as the model constant for turbulence modeling in the larynx.

Rozema (2015) has also shown that after a Taylor expansion of τij in Eq. (5), consistency between AMD and the exact subgrid-scale stress tensor, τij, can be proved. After Taylor expansion of the eddy dissipation of the exact subgrid-scale stress tensor, τ i j S i j, consistency with R5 in Eq. (17) can be also proved. If the exact eddy dissipation gives zero dissipation, then the term R5 gives zero dissipation as well. This means that the AMD model can be switched off for flows where the exact eddy dissipation is vanishing. Thus, the AMD model also switches off when no subgrid energy is created (Rozema , 2015; Vreugdenhil and Taylor, 2018).

The geometry of vocal folds is based on the M5 parametric shape by Scherer (2001); see Fig. 1. The false vocal folds were specified according to data published by Agarwal (2003). The geometrical model is 3D, having a square cross section at inlet 12 × 12 mm. Details can be found in Šidlof (2015).

FIG. 1.

(Color online) Mid-coronal (x-y) section of the CFD domain. Divergent (and also starting) position of vocal folds during phonation is shown. The z-normal (front and back) boundaries belong to Γ wall.

FIG. 1.

(Color online) Mid-coronal (x-y) section of the CFD domain. Divergent (and also starting) position of vocal folds during phonation is shown. The z-normal (front and back) boundaries belong to Γ wall.

Close modal

The boundary conditions for the CFD model are summarized in Table I. The flow is driven by constant pressure difference, P k = p ¯ / ρ = 300 m 2 / s s, between the inlet Γ in and outlet Γ out. The velocity on Γ in and Γ out is computed from the flux. Turbulence initialization at inlet is not used because turbulent intensity upstream of the glottis is low and the flow at inlet can be considered laminar (Lasota and Šidlof, 2019).

TABLE I.

Boundary conditions for the filtered flow velocity, u ¯, and static pressure, p ¯. The symbol n is the unit outer normal and h ( x , t ) is the prescribed displacement of the vocal folds.

u ¯ ( m / s ) p ¯ ( Pa )
Γ in  u ¯ = 0  if   u ¯ n < 0 ,  350 
  u ¯ n = 0  if   u ¯ n > 0    
Γ out  u ¯ = 0  if   u ¯ n < 0 , 
   u ¯ n = 0  if   u ¯ n > 0    
Γ bVF , Γ uVF  u ¯ 1 = 0 , u ¯ 2 = t h ( x , t ) ,  p ¯ n = 0 
  u ¯ 3 = 0   
Γ wall  u ¯ = 0  p ¯ n = 0 
u ¯ ( m / s ) p ¯ ( Pa )
Γ in  u ¯ = 0  if   u ¯ n < 0 ,  350 
  u ¯ n = 0  if   u ¯ n > 0    
Γ out  u ¯ = 0  if   u ¯ n < 0 , 
   u ¯ n = 0  if   u ¯ n > 0    
Γ bVF , Γ uVF  u ¯ 1 = 0 , u ¯ 2 = t h ( x , t ) ,  p ¯ n = 0 
  u ¯ 3 = 0   
Γ wall  u ¯ = 0  p ¯ n = 0 

On the moving boundaries, Γ bVF and Γ uVF, the flow velocity is equal to the velocity of the moving vocal fold surface, given by the function h ( x , t ). The function h ( x , t ) based on the sinusoidal displacement, w 1 , 2 = A 1 , 2 sin ( 2 π f o t + ξ 1 , 2 ), ensures the vibrating motion of vocal folds in the medial-lateral (y) direction with two degrees of freedom. In the current simulation, the vocal folds oscillate symmetrically with a frequency fo = 100 Hz, amplitudes at the superior and inferior vocal fold margin are A 1 , 2 = 0.3 mm, and the phase difference is ξ 1 ξ 2 = π / 2 between the inferior and superior vocal fold margin. During vocal fold oscillation cycle, the medial surface convergence angle, ψ / 2 (see Fig. 1), alternates between ± 10 °. In this study, the oscillation of the vocal folds allows closing/opening the glottal gap, g, in the range 0.42–1.46 mm.

In wall-bounded flows, the presence of solid walls fundamentally influences the flow dynamics, turbulence generation, and transport in the near-wall regions due to significant viscous stresses. The accuracy of the numerical simulation is, thus, closely related to the grid resolution near the fixed walls. According to the classification by Pope (2000), LESs of wall-bounded flows can be classified as large-eddy simulations with near-wall resolution (LES-NWR) with a grid sufficiently fine to resolve 80% of the turbulent energy in the boundary layer, and large-eddy simulation with near-wall modeling (LES-NWM), which employs a modeling approach similar to Reynolds-averaged Navier-Stokes (RANS) in the near-wall region. For these simulations, an important parameter is the wall unit, y + = u τ y / ν, where u τ is the friction velocity and y is the dimensional distance in normal direction from the wall. Using the same normalization, x+ and z+ denote the dimensionless streamwise and spanwise distances, respectively. Wall units are also commonly used to indicate LES adequacy. According to Georgiadis (2010) and Jiang and Lai (2016), in LES-NWR, the theoretical limits for the grid spacing adjacent to the wall are 50 x + 150 , y + < 1 and 15 z + 40, with at least 3–5 gridpoints for 0 < y + < 10. The computational mesh in the current CFD simulation is block structured to capture well the boundary layer and consists of 2.2 M hexahedral elements (Lasota , 2021). An open-source 3D finite volume mesh generator, blockMesh (part of OpenFOAM), was used to build the mesh. The mesh deforms in time due to vocal fold oscillation. On the boundary, Γ bVF, at the critical time when the vocal folds are maximally adducted were evaluated values y avg + = 1.77 , z + = 14 and x + = 8.

The Navier-Stokes equations were discretized using the collocated cell-centered finite volume method. Fletcher (2012) demonstrated that even-ordered derivatives in the truncation error are associated with numerical dissipation, and odd-ordered spatial derivatives are associated with the numerical dispersion in the solution. Ideally, LES simulations should use schemes with low numerical dissipation. The non-dissipative central differencing scheme, which was applied in this study, allows an accurate representation of the changing flow field (Launchbury, 2016). The discretization of the diffusion term is split into an orthogonal and cross-diffusion term using a procedure described in Jasak (1996). Unlike the discretization of the temporal, convective, and orthogonal parts of the diffusive term, the non-orthogonal correctors are treated explicitly.

CFD simulations were run in parallel on the computational cluster Charon [Metacentrum NGI (National Grid Infrastructure)—Technical University of Liberec, using 20 Intel Xeon Silver 4114 cores; Santa Clara, CA] and SGI (silicon graphics based on Intel processors) Altix UV (ultraviolet) 100 supercomputer Fox (Computing center of the Czech Technical University in Prague, using 20 Intel Xeon Nehalem cores).

To have sufficient resolution in the spectrum of the aeroacoustic signal, a sufficiently long simulation time of 20 periods of vocal fold vibration is needed ( t = 0.2 s). For these settings, one CFD simulation required 27–37 days (about 15 000 core-hours of computational time).

In this part, CFD simulations using different turbulence modeling approaches are presented; see Table II. The CFD model validation performed on the periodic hill and backward-facing step is reported in Lasota (2022).

TABLE II.

Overview of the CFD simulations.

Case Type SGS model Cluster Wall-time
LAM  Laminar  —  Charon  27 days 13 h 
OE  LES  OE  Charon  34 days 5 h 
WALE  LES  WALE  Charon  37 days 13 h 
AMD  LES  AMD  Fox  34 days 18 h 
Case Type SGS model Cluster Wall-time
LAM  Laminar  —  Charon  27 days 13 h 
OE  LES  OE  Charon  34 days 5 h 
WALE  LES  WALE  Charon  37 days 13 h 
AMD  LES  AMD  Fox  34 days 18 h 

1. Laryngeal flow rate

Figure 2 shows the glottal opening and flow rates during the last four simulated cycles of vocal fold oscillation. The time, tN, corresponds to the instant when the inferior margins of the vocal folds approach most and reduce the glottal opening to 5.58 mm 2 ( g = 0.465 mm). Time instant, tC, is the maximum approach of the superior vocal fold margins, where the glottal opening drops to 4.98 mm 2 ( g = 0.415 mm). The third time instant, tO, corresponds to the maximum glottal opening of 17.51 mm 2 ( g = 1.459 mm).

FIG. 2.

(Color online) Glottal area waveforms (GAW) and flow rates, Q, during four oscillation cycles. Time instants for further analysis, t N = 0.1900 s , t C = 0.1927 s, and t O = 0.1963 s.

FIG. 2.

(Color online) Glottal area waveforms (GAW) and flow rates, Q, during four oscillation cycles. Time instants for further analysis, t N = 0.1900 s , t C = 0.1927 s, and t O = 0.1963 s.

Close modal

The subgrid-scale models affected the flow rates Q [l/s] (see Fig. 2): the predicted peak flow rate in the laminar case is higher than that in the OE, WALE, and AMD subgrid-scale models by 16.76%, 5.26%, and 9.3%, respectively. This is caused by the different values of the turbulent viscosity, νt, which adds to the molecular viscosity and limits the flow rate through the glottal constriction. The laminar model does not capture the influence of small-scale turbulence, which corresponds to ν t = 0. The WALE and OE compute with nonzero turbulent viscosity, with the latter significantly higher due to the already mentioned deficiency of the OE, which overestimates the turbulent viscosity near the vocal fold surfaces. The minimum and maximum flow rates are 0.122 l / s by using OE and 0.358–0.434 l/s by using the WALE and LAM models, respectively. The peak flow rate predicted by the AMD model occurs slightly sooner than in other simulations.

2. Vorticity field

Vorticity ( ω = × u ¯) is commonly used for characterizing turbulent structures in cases with no entrainment rotation. The vorticity fields reveal the shear layers, where vortices are shed as a consequence of Kelvin-Helmholtz instability. The vortices may undergo successive instabilities, leading to a direct kinetic-energy cascade toward the small scales.

Figure 3 shows vorticity fields presented in the mid-coronal plane (x-y). The supraglottal jet deflects stochastically toward either of the ventricular folds. This behavior is not a consequence of the subgrid-scale model, it is caused by the bistability of the flow in this symmetric geometry (Erath and Plesniak, 2010; Lodermeyer , 2015). Detailed analysis of the vorticity within the glottal region shows that the average value of vorticity in the glottal region is similar for all of the subgrid-scale models.

FIG. 3.

(Color online) Vorticity fields, | ω |, in the mid-coronal plane in range (0,30 000) (1/s).

FIG. 3.

(Color online) Vorticity fields, | ω |, in the mid-coronal plane in range (0,30 000) (1/s).

Close modal

Figure 4 shows a complementary view on the magnitude of the vorticity vector, | ω |, in the midsagittal plane (x-z). The simulation with the AMD model predicts low vorticity near the glottis. The absence of vorticity may imitate the situation in the realistic larynx where the jet is frequently stopped and renewed and, thus, the turbulent eddies are forced to be dissipated.

FIG. 4.

(Color online) Vorticity fields, | ω |, in the midsagittal plane in range (0,30 000) (1/s).

FIG. 4.

(Color online) Vorticity fields, | ω |, in the midsagittal plane in range (0,30 000) (1/s).

Close modal

3. Turbulent viscosity field

The effect of the unresolved turbulent subgrid scales on the resolved scales is carried by the turbulent viscosity, νt, represented by Eq. (8) for ν t O, Eq. (9) for ν t W and Eq. (17) for ν t A.

Figure 5 shows that the turbulent viscosity predicted by the simulation with the OE model is very high in regions of pure shear, especially within the glottis. This is probably the reason why the simulation with the OE model usually predicted the lowest intraglottal velocity. In contrast to this, WALE and AMD subgrid-scale models predicted considerably lower turbulent viscosity in the shear layers at t C. The fields computed by the AMD model seem to be similar to fields computed by WALE with spots of gently higher turbulent viscosity at t C. The other situation occurs at t O when the turbulent viscosity predicted by the AMD model is around 2× higher than that predicted by OE and 5× higher than that predicted by WALE.

FIG. 5.

(Color online) Turbulent viscosity, νt ( m 2 s 1 ), in the mid-coronal plane at tC and tO.

FIG. 5.

(Color online) Turbulent viscosity, νt ( m 2 s 1 ), in the mid-coronal plane at tC and tO.

Close modal

Figure 6 shows turbulent viscosity fields in the midsagittal plane. The simulation with OE predicted 2× higher turbulent viscosity located at the vicinity of the inferior margin of the vocal folds than that at others. The narrow barrier of turbulent viscosity at tN in the case with AMD reduced the flow rate just by 0.8% (1 ml/s) compared to the case with WALE.

FIG. 6.

(Color online) Turbulent viscosity, νt ( m 2 s 1 ), in the midsagittal plane at tC and tO.

FIG. 6.

(Color online) Turbulent viscosity, νt ( m 2 s 1 ), in the midsagittal plane at tC and tO.

Close modal
The CAA model is employed after the CFD simulation is completed and the results are stored on the hard disk. The CAA model is based on the scalar perturbed convective wave equation for the acoustic potential, ψa (Hüppe , 2014),
1 c 0 2 D 2 ψ a D t 2 ( ψ a ) = 1 ρ 0 c 0 2 D p i c D t ,
(19)
where D / D t = / t + u 0 is the material derivative. The overbars ( ¯) from the LES filtering were dropped for simplicity. The speed of sound, c 0 = 351.31 m / s, and ambient density, ρ 0 = 1.1493 kg / m 3, correspond to the temperature of exhaled air of 34 ° C (Anghel and Iacobescu, 2013). The right-hand-side term of Eq. (19) is the aeroacoustic sound source, computed from the incompressible pressure, pic, obtained from the CFD simulation. Finally, the acoustic pressure can be calculated using the relation p a = ρ 0 D ψ a / D t. In our case, there is no background flow, u 0, in front of the lips and within the upper part of the vocal tract. Therefore, p a = ρ 0 ψ a / t can be used because the convective part can be neglected. More details on the CAA model can be found in Lasota (2022).

Figure 7 shows the geometry used for aeroacoustic simulations of vowel /ɑ/, where they are (from left to right) the perfectly matched layer (PML) with a thickness of 0.3 cm at the inlet, larynx with vocal folds, and false vocal folds (2 cm), vocal tract (vowels /ɑ, o/, 17.46 cm; /u/, 18.25 cm; /i, æ/, 16.67 cm) and the radiation zone (RZ) protected by 2–3 PMLs with a thickness of 0.5 cm. PML is a technique published first by Berenger (1994); the original method was modified by Kaltenbacher (2013) by adding damping layers to guarantee that no wave reflections occur at boundaries.

FIG. 7.

(Color online) Acoustic mesh composed of the perfectly matched layers (PML), larynx, vocal tract, and radiation zone (RZ).

FIG. 7.

(Color online) Acoustic mesh composed of the perfectly matched layers (PML), larynx, vocal tract, and radiation zone (RZ).

Close modal

The geometry of the vocal tract was modeled from frustums (0.397 cm long) that were concatenated one after another. The shape of the frustums was defined according to the vocal tract area function measured by MRI (Story , 1996). About 33 000 hexahedral first-order finite elements were used for each vowel.

The vocal tract was conformally attached to the larynx. The connection is formed by two layers of hexahedral cells with minor influence on wave propagation. The right patch of the vocal tract was attached to the RZ.

The element length, Δ l a, of the coarse acoustic mesh and time step, Δ t a, for the aeroacoustic simulation are given by estimations, Δ l a c 0 / ( 20 f max ) and Δ t a 1 / ( 20 f max ), assuming that 20 linear finite elements per 1 acoustic wavelength are sufficient (Hüppe, 2012). In this case, the spatial discretization is limited by 3.43 mm and time step is limited by 1 × 10 5 s to resolve properly acoustic frequencies up to f max = 5 kHz.

The partial differential equation (19) for the acoustic potential, ψa, which is solved numerically in the acoustic domain, is equipped with zero initial conditions and boundary condition ψ a n = 0 , where n is the outward unit normal. This boundary condition can be interpreted as a perfect reflection of a sound wave from a barrier; the condition is also called “sound hard” (Schoder, 2018). In following CAA simulations, the sound hard condition is applied at all of the solid boundaries except the inflow and outflow, where PML is used.

The numerical solution of the aeroacoustic problem proceeds in the following steps (see Fig. 8):

  • The unsteady flow field in the larynx is computed in OpenFOAM on a fine CFD mesh over 20 periods of vocal fold oscillation;

  • the sound sources in the larynx are computed by openCFS and conservatively interpolated onto the coarse CAA mesh, whereas the interpolation back to the fine original CFD mesh is used to visualize the results in high resolution (Figs. 9–12), and

  • in the last step, the aeroacoustic sources are conservatively interpolated onto the coarse CAA mesh and the wave propagation is simulated by openCFS.

FIG. 8.

(Color online) Visualization of data interpolation during the whole workflow.

FIG. 8.

(Color online) Visualization of data interpolation during the whole workflow.

Close modal
FIG. 9.

(Color online) Sound sources, D p i c / D t, from Eq. (19) at time instant, tC. Twenty iso-surfaces in the range ± 2 × 10 5 Pa / s are displayed as 3D (positive sources, purple; negative, green).

FIG. 9.

(Color online) Sound sources, D p i c / D t, from Eq. (19) at time instant, tC. Twenty iso-surfaces in the range ± 2 × 10 5 Pa / s are displayed as 3D (positive sources, purple; negative, green).

Close modal
FIG. 10.

(Color online) Sound sources, D p i c / D t, from Eq. (19) at time instant. tO, but mapped to the mesh with closed-divergent vocal folds positions (tC). Twenty iso-surfaces in the range ± 2 × 10 5 Pa / s are displayed as 3D (positive sources, purple; negative, green).

FIG. 10.

(Color online) Sound sources, D p i c / D t, from Eq. (19) at time instant. tO, but mapped to the mesh with closed-divergent vocal folds positions (tC). Twenty iso-surfaces in the range ± 2 × 10 5 Pa / s are displayed as 3D (positive sources, purple; negative, green).

Close modal
FIG. 11.

(Color online) Sound sources in the mid-coronal plane at the fundamental frequency, f o = 100 Hz, and harmonic frequency, f 9 = 1000 Hz.

FIG. 11.

(Color online) Sound sources in the mid-coronal plane at the fundamental frequency, f o = 100 Hz, and harmonic frequency, f 9 = 1000 Hz.

Close modal
FIG. 12.

(Color online) Sound sources at a nonharmonic frequency, f = 1235 Hz.

FIG. 12.

(Color online) Sound sources at a nonharmonic frequency, f = 1235 Hz.

Close modal

The computational time needed for one CAA simulation is much lower than that needed for one CFD simulation, about 5 h on a single central processing unit (CPU) core compared to 27–37 days on 20 cores. The conservative interpolation of sound sources (Schoder , 2019; Schoder , 2020; Schoder , 2021b) from the CFD mesh to the CAA mesh was performed by the cfsdat tool, which is part of the openCFS library.

1. Sound sources (time domain)

The distribution of the aeroacoustic sources in the computational domain covering the larynx varies throughout the vocal fold oscillation period. The jet is surrounded by spots of strong positive and negative acoustic sources related to turbulent eddies created from shear layers of the jet; see also Schoder (2021a). Figure 9 shows the aeroacoustic sound source distribution corresponding to the closed-convergent position of vocal folds, which are visualized on the CFD mesh with the closed-divergent position of vocal folds. The aeroacoustic sources computed on the moving geometry with oscillating vocal folds are mapped to a fixed geometry. This is performed because the current version of the acoustic solver cannot handle moving meshes.

Figure 10 shows sound source distribution when vocal folds are fully open, and the results are again mapped to the closed-divergent domain. Strong sound sources within the glottis and 1 cm from the glottis are observed. The 3D view shows uniform distribution of sound sources in the glottis. This can be caused by the presence of high turbulent viscosity.

2. Sound sources (frequency domain)

The conversion from the time to frequency domain was made by the field fast Fourier transform (field FFT), which brings a better insight into the spatial distribution of the sound sources at distinct frequencies. Figure 11 shows that the results obtained from the LES without subgrid-scale model (LAM) show slightly more intensive sound sources within the glottis at f o = 100 Hz than that at other cases. This correlates with the flow rate amplitude, which is also higher in the LAM case (see Fig. 2). On the other side, the intensity of sound sources at f = 1000 Hz is opposite, and AMD is 2.5–4× higher compared to the others. At these higher frequencies, dominant sound sources do not occur within the glottis but in the places where the fast glottal jet interacts with the ventricular folds and the slowly moving recirculating air in the supraglottal volume.

Figure 12 shows sound sources at a random nonharmonic frequency, f = 1235 Hz, where the sound sources are distributed further downstream of the glottis. The integrity of sound sources can be observed along with several local sound spots at the superior edge of the vocal folds.

3. Wave propagation (time domain)

In the following, the acoustic pressure fields, p a ( x , t ), will be analyzed as a solution of Eq. (19). The acoustic pressure has been tracked at distances of 1 and 16 cm from the lips to create audio recordings; see Fig. 13. Two PML layers around the RZ damped the acoustic waves on walls properly as no reflections are observed.

FIG. 13.

(Color online) Acoustic wave propagation in the radiation field after three periods of vocal folds oscillation.

FIG. 13.

(Color online) Acoustic wave propagation in the radiation field after three periods of vocal folds oscillation.

Close modal

Table III compares individual cases in terms of p rms a = ( 1 / T ) 0 T ( p a ) 2 d t and SPL = 20 log 10 ( p rms a / p ref a ), where p ref a = 20 μ Pa is the hearing threshold. The total sound pressure level (SPL) value of the radiated sound can be used as a measure of the acoustic energy transferred from the larynx through the vocal tract. Minor differences in SPL are observed between the back close vowel, /u/, and front open vowel, /ɑ/. The simulations of front open/mid vowel, /æ/, transferred most energy of all of the simulated vowels. The simulations based on the AMD model predicted the highest or similar SPLs of all of the vowels, except the front-close vowel, /i/. This may lead to the conclusion that WALE model can transfer more energy in simulations with the close vowels.

TABLE III.

Root-mean-square acoustic pressures, p rms a (Pa), and total sound pressure levels, SPLs (dB), tracked at 16 cm.

Case p rms a SPL p rms a SPL p rms a SPL
LAM  ɑ  0.0031  43.69  æ  0.0052  48.37  0.0026  42.12 
OE  ɑ  0.0016  37.89  æ  0.0021  40.45  0.0011  35.16 
WALE  ɑ  0.0028  42.80  æ  0.0049  47.84  0.0022  40.68 
AMD  ɑ  0.0032  44.13  æ  0.0051  48.12  0.0016  38.01 
LAM  0.0019  39.48  0.0019  39.69       
OE  0.0016  38.31  0.0013  36.58       
WALE  0.0018  38.89  0.0017  38.43       
AMD  0.0026  42.23  0.0018  39.22       
Case p rms a SPL p rms a SPL p rms a SPL
LAM  ɑ  0.0031  43.69  æ  0.0052  48.37  0.0026  42.12 
OE  ɑ  0.0016  37.89  æ  0.0021  40.45  0.0011  35.16 
WALE  ɑ  0.0028  42.80  æ  0.0049  47.84  0.0022  40.68 
AMD  ɑ  0.0032  44.13  æ  0.0051  48.12  0.0016  38.01 
LAM  0.0019  39.48  0.0019  39.69       
OE  0.0016  38.31  0.0013  36.58       
WALE  0.0018  38.89  0.0017  38.43       
AMD  0.0026  42.23  0.0018  39.22       

4. Wave propagation (frequency domain)

This section deals with frequency spectra of simulated signals. This kind of analysis can highlight fundamental, harmonic, and nonharmonic frequencies, which are accompanied by broadband noise. The FFT analyses were performed on the signal spanning over 20 periods of the vocal fold oscillation (1 period = 10 ms = 1000 samples) and, thus, the frequency resolution is Δ f = 5 Hz.

Influence of the end of the vocal tract (lips) was tested first, comparing two signals at 1 and 16 cm; see Fig. 14. At 16 cm, the SPLs on all of the frequency components are uniformly lower by 19 dB. Our model of phonation prescribing 350 Pa at inlet of the larynx can be considered to be quiet phonation while, in the other case, it is recommended to use a probe for tracking a signal at least 3 cm from the lips. The envelope overlaying black peaks corresponds to the linear predictive coding (LPC) curve (Pavlidi , 2013), which is widely used to provide a smooth spectral envelope of audio and voice recordings for the purpose of finding positions of formants.

FIG. 14.

(Color online) Two frequency spectra of vocalization of /ɑ/ at the distance of 1 cm and 16 cm from the lips.

FIG. 14.

(Color online) Two frequency spectra of vocalization of /ɑ/ at the distance of 1 cm and 16 cm from the lips.

Close modal

In the following, the frequency spectra at the distance of 16 cm will be analyzed vowel by vowel, and the influence of the subgrid-scale models will be investigated. Note that the LES without subgrid-scale model (LAM) should be considered as the least accurate case, which disregards the effect of small-scale turbulence. However, this approach is often still used in the modeling of laryngeal flow, and the comparison against more accurate LESs is interesting.

a. Vowel /ɑ/.

The spectrum is presented in Fig. 15, where SPLs at fundamental frequency, f o = 100 Hz, and higher harmonics, f 1 = 200 Hz , f 2 = 300 Hz, and so forth are well visible. Peaks at fo are lower than those at f1 and f2. If the signal is processed immediately behind the vocal folds, the SPL at fo is definitely the highest, whereas the signal tracked in RZ includes the effect of the vocal tract for vowel /ɑ/. The peaks (formants) are key parameters, which define a vowel. The close distance between the first two formants, F 1 F 2, is typical for vowels /ɑ, o, u/. It should be noted that the first two formants are important for vowel distinctness. The AMD model resulted in the higher peak at the first formant around 900 Hz by 8 dB and, thus, contributes to the better vowel clarity compared to other models.

FIG. 15.

(Color online) Four frequency spectra of vocalization of /ɑ/ at 16 cm from the lips.

FIG. 15.

(Color online) Four frequency spectra of vocalization of /ɑ/ at 16 cm from the lips.

Close modal
b. Vowel /æ/.

Figure 16 shows the frequency spectrum for this front and open vowel, /æ/, and corresponds to the highest total SPL from all of the simulated vowels, up to 48 dB.

FIG. 16.

(Color online) Four frequency spectra of vocalization of /æ/ at 16 cm from the lips.

FIG. 16.

(Color online) Four frequency spectra of vocalization of /æ/ at 16 cm from the lips.

Close modal
c. Vowel /i/.

The frequency spectrum for this front and close vowel is plotted in Fig. 17. The spectrum shows that the simulations with OE predicted the lowest total SPL (around 35 dB), which is related to the well-known overpredicting of turbulent viscosity by the OE model and can be observed in all of the spectra.

FIG. 17.

(Color online) Four frequency spectra of vocalization of /i/ at 16 cm from the lips.

FIG. 17.

(Color online) Four frequency spectra of vocalization of /i/ at 16 cm from the lips.

Close modal
d. Vowel /o/.

Figure 18 shows the frequency spectrum for this back and mid-close vowel for which it is typical that F 2 F 3 are far apart. Simulations based on WALE predicted the most visible third formant of all of the vowels. To distinguish vowels, human ears are most sensitive to the two lowest formants. For that reason, the AMD model can be more valuable for voice research.

FIG. 18.

(Color online) Four frequency spectra of vocalization of /o/ at 16 cm from the lips.

FIG. 18.

(Color online) Four frequency spectra of vocalization of /o/ at 16 cm from the lips.

Close modal
e. Vowel /u/.

The frequency spectrum for this back and close vowel is in Fig. 19. The spectrum shows that the second formant of the LAM case at 1100 Hz is higher by 8 dB compared to the second formants of the WALE and AMD. However, the audio recording from the simulation with LAM is subjectively not much clearer compared to the others.

FIG. 19.

(Color online) Four frequency spectra of vocalization of /u/ at 16 cm from the lips.

FIG. 19.

(Color online) Four frequency spectra of vocalization of /u/ at 16 cm from the lips.

Close modal

The formant locations, F 1 F 2, computed from local maxima of the LPC curves in Figs. 15–19, are compared with measurements of natural speech by Ireland (2015), Peterson and Barney (1952), and Story (1996) in Fig. 20. All of the simulated vowels lie inside of the measured ranges, which confirms the usability of the acoustic grids based on the circular cross sections.

FIG. 20.

(Color online) Formant ranges measured by Peterson and Barney (1952) and formant locations obtained by Story (1996) and Ireland (2015) are compared with the simulated formants.

FIG. 20.

(Color online) Formant ranges measured by Peterson and Barney (1952) and formant locations obtained by Story (1996) and Ireland (2015) are compared with the simulated formants.

Close modal

Readers can assess the clarity of simulated vowels by downloading the audio recordings.2 Please note that for comfortable, listening recordings at the distance of 1 cm from the lips were used.

In this article, models of phonation based on LESs with subgrid-scale models contributing to the solution of turbulent flow were presented. It was shown that the OE model overpredicts the turbulent viscosity in the boundary layer adjacent to the vocal folds and shear layers of the glottal jet (Fig. 5). The difference in glottal flow rate among the simulations is clearly induced by the subgrid-scale model, which adds the turbulent viscosity to the molecular viscosity of air and hinders the airflow in the glottis (Fig. 2). The WALE model produced zero eddy viscosity in cases of pure shear flow (Fig. 6) and, hence, the flow simulation with WALE predicted by 5% higher the maximum transglottal flow rate than did AMD. Despite this fact, the phonation simulation based on the AMD model transferred more or equal energy in terms of total SPL than WALE for all of the vowels except the front-close vowel, /i/ (Table III). The WALE model, which is known to handle turbulent viscosity at the near-wall and high-shear regions more precisely than the OE model (Figs. 5 and 6), resulted in higher total SPLs than OE in all of the cases (Table III). The OE model gives acceptable results, in general, but peaks of frequency formants are hardly visible and weaker compared to the WALE or AMD models (Figs. 16–19). The WALE model resulted in higher third formants in high-frequency bandwidth, most of all, for the subgrid-scale models (Figs. 16–19). However, the third formant is not crucial for vowel characterization. On the other side, the AMD model resulted in slightly higher SPLs at harmonic and formant frequencies up to the second formant for all of the vowels (Figs. 16–19). First formants are higher in cases /ɑ, æ/ at least by 8 dB (Figs. 15 and 16). These positive findings can be attributed to beneficial features of the AMD model: consistency with the exact subgrid-scale stress tensor, τij, no requirements on the approximation of the LES filter width, and usability on anisotropic meshes.

Based on recordings of simulated vowels, it can be concluded that the model of phonation employing the AMD subgrid-scale model is much closer to natural speech than large-eddy simulations with conventional subgrid-scale models.

The research was supported by the Student Grant Scheme at the Technical University of Liberec through Project No. SGS-2022-3016. The Graz group acknowledges support from the ÖAW research grant “Understanding voice disorders,” received from Dr. Anton Oelzelt-Newin'sche Stiftung.

1

See https://gitlab.com/mlasota/myFoam (Last viewed January 21, 2023).

2

See supplementary material at https://www.scitation.org/doi/suppl/10.1121/10.0017202 to download the (.wav) recordings of simulated vowels sorted by different subgrid-scale models.

1.
Agarwal
,
M.
,
Scherer
,
R.
, and
Hollien
,
H.
(
2003
). “
The false vocal folds: Shape and size in frontal view during phonation based on laminagraphic tracings
,”
J. Voice
17
(
2
),
97
113
.
2.
Anghel
,
M.-A.
, and
Iacobescu
,
F.
(
2013
). “
The influence of temperature and co2 in exhaled breath
,” in
16th International Congress of Metrology
,
EDP Sciences
, p.
10012
.
3.
Avhad
,
A.
,
Li
,
Z.
,
Wilson
,
A.
,
Sayce
,
L.
,
Chang
,
S.
,
Rousseau
,
B.
, and
Luo
,
H.
(
2022
). “
Subject-specific computational fluid-structure interaction modeling of rabbit vocal fold vibration
,”
Fluids
7
(
3
),
97
.
4.
Berenger
,
J.-P.
(
1994
). “
A perfectly matched layer for the absorption of electromagnetic waves
,”
J. Comput. Phys.
114
(
2
),
185
200
.
5.
Bodaghi
,
D.
,
Jiang
,
W.
,
Xue
,
Q.
, and
Zheng
,
X.
(
2021
). “
Effect of supraglottal acoustics on fluid-structure interaction during human voice production
,”
J. Biomech. Eng.
143
(
4
),
041010
.
6.
Döllinger
,
M.
,
Kobler
,
J. A.
,
Berry
,
D.
,
Mehta
,
D.
,
Luegmair
,
G.
, and
Bohr
,
C.
(
2011
). “
Experiments on analysing voice production: Excised (human, animal) and in vivo (animal) approaches
,”
Curr. Bioinform.
6
(
3
),
286
304
.
7.
Erath
,
B.
, and
Plesniak
,
M.
(
2010
). “
An investigation of asymmetric flow features in a scaled-up driven model of the human vocal folds
,”
Exp. Fluids
49
(
1
),
131
146
.
8.
Falk
,
S.
,
Kniesburges
,
S.
,
Schoder
,
S.
,
Jakubaß
,
B.
,
Maurerlehner
,
P.
,
Echternach
,
M.
,
Kaltenbacher
,
M.
, and
Döllinger
,
M.
(
2021
). “
3D-FV-FE aeroacoustic larynx model for investigation of functional based voice disorders
,”
Front. Physiol.
12
,
226
.
9.
Ferziger
,
H.
(
1998
). “
Direct and large eddy simulation of turbulence
,”
Numer. Methods Fluid Mech.
16
,
53
73
.
10.
Ffowcs Williams
,
J. E.
, and
Hawkings
,
D. L.
(
1969
). “
Sound generation by turbulence and surfaces in arbitrary motion
,”
Philos. Trans. R. Soc. London, Ser. A, Math. Phys. Sci.
264
(
1151
),
321
342
.
11.
Fletcher
,
C. A. J.
(
2012
). Computational Techniques for Fluid Dynamics 2: Specific Techniques for Different Flow Categories (Springer Science & Business Media).
12.
Georgiadis
,
N. J.
,
Rizzetta
,
D. P.
, and
Fureby
,
C.
(
2010
). “
Large-eddy simulation: Current capabilities, recommended practices, and future research
,”
AIAA J.
48
(
8
),
1772
1784
.
13.
Hüppe
,
A.
(
2012
). “
Spectral finite elements for acoustic field computation
,” Ph.D. thesis,
Alps-Adriatic University of Klagenfurt
.
14.
Hüppe
,
A.
,
Grabinger
,
J.
,
Kaltenbacher
,
M.
,
Reppenhagen
,
A.
,
Dutzler
,
G.
, and
Kühnel
,
W.
(
2014
). “
A non-conforming finite element method for computational aeroacoustics in rotating systems,” in
20th AIAA/CEAS Aeroacoustics Conference
, p.
2739
.
15.
Ireland
,
D.
,
Knuepffer
,
C.
, and
McBride
,
S. J.
(
2015
). “
Adaptive multi-rate compression effects on vowel analysis
,”
Front. Bioeng. Biotechnol.
3
,
118
.
16.
Jasak
,
H.
(
1996
). “
Error analysis and estimation for the finite volume method with applications to fluid flows
,” Ph.D. thesis,
Imperial College London
.
17.
Jiang
,
X.
, and
Lai
,
C.-H.
(
2016
).
Numerical Techniques for Direct and Large-Eddy Simulations
(
CRC Press
,
Boca Raton, FL
).
18.
Kaltenbacher
,
B.
,
Kaltenbacher
,
M.
, and
Sim
,
I.
(
2013
). “
A modified and stable version of a perfectly matched layer technique for the 3-d second order wave equation in time domain with an application to aeroacoustics
,”
J. Comput. Phys.
235
,
407
422
.
19.
Kniesburges
,
S. L.
,
Thomson
,
S.
,
Barney
,
A.
,
Triep
,
M.
,
Šidlof
,
P.
,
Horáček
,
J.
,
Brücker
,
C.
, and
Becker
,
S.
(
2011
). “
In vitro experimental investigation of voice production
,”
Curr. Bioinform.
6
(
3
),
305
322
.
20.
Lasota
,
M.
(
2022
). “
Large-eddy simulation for aeroacoustics of human phonation
,” Ph.D. thesis,
TU Liberec
.
21.
Lasota
,
M.
, and
Šidlof
,
P.
(
2019
). “
Large-eddy simulation of flow through human larynx with a turbulence grid at inlet,” in
Topical Problems Fluid Mechanics
.
22.
Lasota
,
M.
,
Šidlof
,
P.
,
Kaltenbacher
,
M.
, and
Schoder
,
S.
(
2021
). “
Impact of the sub-grid scale turbulence model in aeroacoustic simulation of human voice
,”
Appl. Sci.
11
(
4
),
1970
.
23.
Launchbury
,
D. R.
(
2016
).
Unsteady Turbulent Flow Modelling and Applications
(
Springer
,
Berlin
).
24.
Lesieur
,
M.
,
Métais
,
O.
, and
Comte
,
P.
(
2005
).
Large-Eddy Simulations of Turbulence
(
Cambridge University Press
,
Cambridge, UK
).
25.
Lodermeyer
,
A.
,
Becker
,
S.
,
Döllinger
,
M.
, and
Kniesburges
,
S.
(
2015
). “
Phase-locked flow field analysis in a synthetic human larynx model
,”
Exp. Fluids
56
(
4
),
1
13
.
26.
Mattheus
,
W.
, and
Brücker
,
C.
(
2011
). “
Asymmetric glottal jet deflection: Differences of two- and three-dimensional models
,”
J. Acoust. Soc. Am.
130
(
6
),
EL373
EL379
.
27.
Mihaescu
,
M.
,
Khosla
,
S. M.
,
Murugappan
,
S.
, and
Gutmark
,
E. J.
(
2010
). “
Unsteady laryngeal airflow simulations of the intra-glottal vortical structures
,”
J. Acoust. Soc. Am.
127
(
1
),
435
444
.
28.
Nicoud
,
F.
, and
Ducros
,
F.
(
1999
). “
Subgrid-scale stress modelling based on the square of the velocity gradient tensor
,”
Flow, Turbul. Combust.
62
(
3
),
183
200
.
29.
Pavlidi
,
D.
,
Griffin
,
A.
,
Puigt
,
M.
, and
Mouchtaris
,
A.
(
2013
). “
Real-time multiple sound source localization and counting using a circular microphone array
,”
IEEE Trans. Audio, Speech, Lang. Process.
21
(
10
),
2193
2206
.
30.
Peterson
,
G. E.
, and
Barney
,
H. L.
(
1952
). “
Control methods used in a study of the vowels
,”
J. Acoust. Soc. Am.
24
(
2
),
175
184
.
31.
Pope
,
S. B.
(
2000
).
Turbulent Flows
(
Cambridge University Press
,
Cambridge, UK
).
32.
Rozema
,
W.
,
Bae
,
H. J.
,
Moin
,
P.
, and
Verstappen
,
R.
(
2015
). “
Minimum-dissipation models for large-eddy simulation
,”
Phys. Fluids
27
(
8
),
085107
.
33.
Sadeghi
,
H.
,
Döllinger
,
M.
,
Kaltenbacher
,
M.
, and
Kniesburges
,
S.
(
2019a
). “
Aerodynamic impact of the ventricular folds in computational larynx models
,”
J. Acoust. Soc. Am.
145
(
4
),
2376
2387
.
34.
Sadeghi
,
H.
,
Kniesburges
,
S.
,
Falk
,
S.
,
Kaltenbacher
,
M.
,
Schützenberger
,
A.
, and
Döllinger
,
M.
(
2019b
). “
Towards a clinically applicable computational larynx model
,”
Appl. Sci.
9
(
11
),
2288
.
35.
Scherer
,
R.
,
Shinwari
,
D.
,
De Witt
,
J.
,
Zhang
,
C.
,
Kucinschi
,
R.
, and
Afjeh
,
A.
(
2001
). “
Intraglottal pressure profiles for a symmetric and oblique glottis with a divergence angle of 10 degrees
,”
J. Acoust. Soc. Am.
109
(
4
),
1616
1630
.
36.
Schoder
,
S.
(
2018
). “
Aeroacoustic analogies based on compressible flow data
,” Ph.D. thesis,
TU Wien
.
37.
Schoder
,
S.
,
Junger
,
C.
,
Weitz
,
M.
, and
Kaltenbacher
,
M.
(
2019
). “
Conservative source term interpolation for hybrid aeroacoustic computations
,” in
25th AIAA/CEAS Aeroacoustics Conference
, p.
2538
.
38.
Schoder
,
S.
,
Maurerlehner
,
P.
,
Wurzinger
,
A.
,
Hauser
,
A.
,
Falk
,
S.
,
Kniesburges
,
S.
,
Döllinger
,
M.
, and
Kaltenbacher
,
M.
(
2021a
). “
Aeroacoustic sound source characterization of the human voice production-perturbed convective wave equation
,”
Appl. Sci.
11
(
6
),
2614
.
39.
Schoder
,
S.
,
Weitz
,
M.
,
Maurerlehner
,
P.
,
Hauser
,
A.
,
Falk
,
S.
,
Kniesburges
,
S.
,
Döllinger
,
M.
, and
Kaltenbacher
,
M.
(
2020
). “
Hybrid aeroacoustic approach for the efficient numerical simulation of human phonation
,”
J. Acoust. Soc. Am.
147
(
2
),
1179
1194
.
40.
Schoder
,
S.
,
Wurzinger
,
A.
,
Junger
,
C.
,
Weitz
,
M.
,
Freidhager
,
C.
,
Roppert
,
K.
, and
Kaltenbacher
,
M.
(
2021b
). “
Application limits of conservative source interpolation methods using a low Mach number hybrid aeroacoustic workflow
,”
J. Theor. Comput. Acoust.
29
(
01
),
2050032
.
41.
Schwarze
,
R.
,
Mattheus
,
W.
,
Klostermann
,
J.
, and
Brücker
,
C.
(
2011
). “
Starting jet flows in a three-dimensional channel with larynx-shaped constriction
,”
Comput. Fluids
48
(
1
),
68
83
.
42.
Šidlof
,
P.
,
Zörner
,
S.
, and
Hüppe
,
A.
(
2015
). “
A hybrid approach to the computational aeroacoustics of human voice production
,”
Biomech. Model. Mechanobiol.
14
(
3
),
473
488
.
43.
Smagorinsky
,
J.
(
1963
). “
General circulation experiments with the primitive equations: I. The basic experiment
,”
Mon. Weather Rev.
91
(
3
),
99
164
.
44.
Story
,
B. H.
,
Titze
,
I. R.
, and
Hoffman
,
E. A.
(
1996
). “
Vocal tract area functions from magnetic resonance imaging
,”
J. Acoust. Soc. Am.
100
(
1
),
537
554
.
45.
Suh
,
J.
, and
Frankel
,
S.
(
2007
). “
Numerical simulation of turbulence transition and sound radiation for flow through a rigid glottal model
,”
J. Acoust. Soc. Am.
121
(
6
),
3728
3739
.
46.
Tokuda
,
I. T.
, and
Shimamura
,
R.
(
2017
). “
Effect of level difference between left and right vocal folds on phonation: Physical experiment and theoretical study
,”
J. Acoust. Soc. Am.
142
(
2
),
482
492
.
47.
Versteeg
,
H. K.
, and
Malalasekera
,
W.
(
2007
).
An Introduction to Computational Fluid Dynamics: The Finite Volume Method
(
Pearson Education
,
London, UK
).
48.
Vreugdenhil
,
C. A.
, and
Taylor
,
J. R.
(
2018
). “
Large-eddy simulations of stratified plane couette flow using the anisotropic minimum-dissipation model
,”
Phys. Fluids
30
(
8
),
085104
.
49.
Yoshinaga
,
T.
,
Van Hirtum
,
A.
, and
Wada
,
S.
(
2017
). “
Multimodal modeling and validation of simplified vocal tract acoustics for sibilant /s/
,”
J. Sound Vib.
411
,
247
259
.
50.
Yoshizawa
,
A.
, and
Horiuti
,
K.
(
1985
). “
A statistically-derived subgrid-scale kinetic energy model for the large-eddy simulation of turbulent flows
,”
J. Phys. Soc. Jpn.
54
(
8
),
2834
2839
.
51.
Zahiri
,
A.-P.
, and
Roohi
,
E.
(
2019
). “
Anisotropic minimum-dissipation (AMD) subgrid-scale model implemented in openfoam: Verification and assessment in single-phase and multi-phase flows
,”
Comput. Fluids
180
,
190
205
.
52.
Zhang
,
Z.
,
Wu
,
L.
,
Gray
,
R.
, and
Chhetri
,
D. K.
(
2020
). “
Three-dimensional vocal fold structural change due to implant insertion in medialization laryngoplasty
,”
PLoS One
15
(
1
),
1
11
.

Supplementary Material