PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
J Fluids Struct. Author manuscript; available in PMC 2010 November 1.
Published in final edited form as:
J Fluids Struct. 2009 November 1; 25(8): 1299–1317.
doi:  10.1016/j.jfluidstructs.2009.08.002
PMCID: PMC2790218
NIHMSID: NIHMS148284

ANALYSIS OF FLOW-STRUCTURE COUPLING IN A MECHANICAL MODEL OF THE VOCAL FOLDS AND THE SUBGLOTTAL SYSTEM

Abstract

An analysis is made of the nonlinear interactions between flow in the subglottal vocal tract and glottis, sound waves in the subglottal system and a mechanical model of the vocal folds. The mean flow through the system is produced by a nominally steady contraction of the lungs, and mechanical experiments frequently involve a ‘lung cavity’ coupled to an experimental subglottal tube of arbitrary or ill-defined effective length L, on the basis that the actual value of L has little or no influence on excitation of the vocal folds. A simple, self-exciting single mass mathematical model of the vocal folds is used to investigate the sound generated within the subglottal domain and the unsteady volume flux from the glottis for experiments where it is required to suppress feedback of sound from the supraglottal vocal tract. In experiments where the assumed absorption of sound within the sponge-like interior of the lungs is small, the influence of changes in L can be very significant: when the subglottal tube behaves as an open-ended resonator (when L is as large as half the acoustic wavelength) there is predicted to be a mild increase in volume flux magnitude and a small change in waveform. However, the strong appearance of second harmonics of the acoustic field is predicted at intermediate lengths, when L is roughly one quarter of the acoustic wavelength. In cases of large lung damping, however, only modest changes in the volume flux are predicted to occur with variations in L.

Keywords: phonation, voiced speech, vocal folds, vortex sound

1. Introduction

Titze (1988) has argued that feedback coupling of sound waves in the upper and lower sections of the vocal tract on the motion of the vocal folds must generally be relatively weak. Such coupling is likely to be at its strongest at higher frequencies when the voicing frequency approaches the first formant (Joliveau et al. 2004). Then voicing sources can interact with resonances to produce undesirable, involuntary and abrupt changes in frequency; this has been observed during phonation, particularly for resonances of the lower (or subglottal) section of the vocal tract (Titze and Story 1995; Austin and Titze 1997; Titze 2008).

Experimental studies of phonation that involve either an excised larynx or a laboratory model of the larynx are frequently designed to avoid unwanted coupling with the subglottal region by use of a long subglottal duct whose resonance frequencies are very much smaller than voiced frequencies of interest (Titze and Story, 1995; Alipour and Scherer, 2001; Thomson et al., 2005; Zhang et al., 2006a). On the other hand, some experiments are performed using physical models that replicate more closely the human subglottal system (e.g. Zhang et al., 2006b, 2009). The length of the subglottal tube upstream of the larynx is then typically between 17 and 20 cm and it usually terminates in a large reservoir or plenum. The whole arrangement is intended to provide a simplified mechanical model of the trachea, bronchi, and lungs (Flanagan, 1958; Ishizaka et al., 1976).

Many reports of these experiments omit completely important details of the laboratory subglottal region, so that likely influences on voicing of subglottal geometry and mechanical properties cannot easily be ascertained, even when these influences could play a prominent role in determining measured output acoustics (Zhang et al., 2006a). Similarly, numerical investigations frequently proceed on the assumption that the subglottal region is a constant pressure source, without resonances. Indeed, these difficulties were recognized by Zhang et al. (2006a) who concluded, from detailed measurements supported by linearized analysis of the coupled motion of the vocal folds and the subglottal sound, that the results of many experiments reported in the literature are actually a consequence of the action on the vocal folds of unreported or unmeasured feedback from subglottal resonances.

Estimates of feedback strength and its influence on glottis oscillations have generally involved ad hoc representations of the coupling between the acoustic field and surface tissue in the glottis (Flanagan and Landgraf 1968; Gupta et al. 1973; Zanartu et al. 2007; Titze 2008). In particular, according to an analysis of Zanartu et al. (2007), using a single-mass model of the glottis, the back reaction of acoustic pressures strongly control the maintenance of self-sustaining oscillations. When the vocal folds are coupled only to the lower vocal tract no self-sustained oscillation is observed for any subglottal configuration, apparently in contradiction to the simple single-mass models suggested by Zhang et al. (2006b). Thus the properties of the subglottal tract can strongly affect phonation, although according to Zanartu et al. (2007) the influence of acoustic feedback from the upper tract is more important. Many of these conclusions may be a consequence of the properties of the particular single-mass model of the vocal folds used in their analysis because, for example, the one-mass model of Howe and McGowan (2009) exhibits self-sustained oscillations even in the absence of acoustic feedback. However, a number of experimental studies (e.g. Austin and Titze (1997), Zhang et al. (2006a, 2006b)) record variations of the fundamental frequency with changes in the subglottal and supraglottal tracts.

In this paper subglottal feedback is examined in detail for a mechanical model of the glottis whose fluid mechanical characteristics can be specified precisely, and where the acoustic and aerodynamic pressures associated with vortex shedding during jet formation are calculated with proper account taken of the local geometry of the moving boundaries. A nonlinear analytical model is developed in terms of the theory of aerodynamic sound (Howe 1998; Howe and McGowan 2007) for a subglottal system often used in experiments. The system consists of a rigid tube of constant cross-section joining a mechanical model of the vocal folds and terminating at its remote end at a reflectionless plenum or duct of much larger cross-section, henceforth referred to as the ‘lungs’. A self-sustaining model, requiring no empirical input, is investigated that makes use of a simple, but aerodynamically consistent single-mass model of the vocal folds (described in detail by Howe and McGowan 2009) in accordance with the notion of nonlinear excitation propounded by Fant (1960, Section A.21). The theory determines (for this particular model) the details of the influence of subglottal acoustics on the unsteady volume flux through the glottis, and the dependence of the wave profile on the length of the subglottal tube, and on conditions at the interface with the lungs, where the absorption of sound within the lung complex is shown to strongly influence the acoustic-vocal fold interaction. The whole motion is driven by a steady contraction of the lung cavity. The investigation complements the linearised study of the same problem by Zhang et al. (2006a).

The analytical problem is formulated and formally solved in Sections 2, 3, with the help of an aeroacoustic Green’s function derived in the Appendix. The calculation of the glottal volume flux is discussed in Section 4, and numerical predictions of the coupled motions of the vocal folds and the subglottal acoustics are then presented (Section 5).

2. Coupling of the subglottal tract and the vocal folds

We consider a simplified, analytical model of the subglottal vocal tract and vocal folds, illustrated overall in Figure 1. In this model the upper tract is absent so that sound generated at the glottis radiates directly into free space. This arrangement is used frequently in phonation experiments (although, of course, the configuration is not a strictly physiological model) and precludes any possible feedback from the upper tract (Zhang et al. 2006a). The subglottal tract (or ‘tube’) is just upstream of the vocal folds and is treated as a uniform, hard walled duct of length [ell] and cross-sectional area An external file that holds a picture, illustration, etc.
Object name is nihms148284ig1.jpg. The duct axis is taken to coincide with the interval −[ell] < x < 0 of the x axis of the rectangular coordinate system x = (x, y, z). The ‘open’ end at x = 0 communicates with the ambient free space via an idealized glottal opening, formed by a simplified ‘single mass’ mechanical model of the vocal folds, which is described below. The tract terminates at its inner end (x = −[ell]) with an abrupt increase in cross-sectional area to An external file that holds a picture, illustration, etc.
Object name is nihms148284ig2.jpg followed by a nominally infinite duct that is taken to model the lung plenum – sound radiating from the subglottal tube into this duct proceeds without reflection to x = −∞, which is taken to be equivalent to absorption by the lossy compliance of the lungs.

Figure 1
Mathematical model of the subglottal vocal region coupled to the glottis.

Voicing in vocal fold experiments is initiated by supplying the plenum with a constant volume flow of air. This mimics the steady contraction of the lung cavity and is represented mathematically by the action of a planar, steady volume source that spans the lower lung plenum at x = −[l with macron] (<−[ell]). The steady flow from this source first enters the subglottal tube as a pressure wave that is subsequently multiply reflected within the tube and modified by vortex sound generated by ‘jetting’ at the glottis. It is required to determine the sound radiated from the glottis (at the free space point x in Figure 1, for example); this is equivalent to determining the unsteady volume flux through the glottis because, from the point of view of an observer in free space, the radiation is necessarily of monopole type. In normal adult speech the voicing frequency f0 ~ 125 Hz is very much smaller than the first formant of the subglottal system (~ 700 Hz), so that resonances of the lower vocal tract would not normally be expected to modify significantly voicing properties.

2.1 Single-mass model of the vocal folds

The volume flux through the glottis is modulated by fluctuations in the glottal cross-sectional area produced by aerodynamic coupling between the flow and the vocal folds. In the original single-mass model (Fant 1960; Flanagan 1972) (Figure 2) the glottal flow was approximated by potential flow through an orifice whose cross-section varied in a prescribed manner with time. Although this potential flow is transformed into a free field jet on emerging from the glottis, irrotational flow fluctuations within the glottis are incapable of supplying the feedback necessary to maintain the folds in oscillatory motion.

Figure 2
Geometry of the single mass model of the vocal folds.

During phonation the glottal cross-section resembles more-or-less an ellipse of large aspect ratio, whose major and minor axes 2a and 2b in Figure 3a vary with time. In this paper this elongated shape is approximated by a rectangle of fixed span 2a, marginally smaller than the diameter of the vocal tract near the glottis, and of time-dependent width 2ζ(t) (Figure 3b). To fix ideas, we suppose that the vocal tract has rectangular cross-section at the glottis, of span [ell]3 ≈ 2a and width 2h, as indicated in Figure 2. The variation of ζ(t) is governed by the aerodynamic interaction of the vocal folds and the unsteady flow through the glottis. This interaction is assumed to occur according to a modified single-mass model introduced by Howe and McGowan (2009); the model includes the influence of quasi-two-dimensional flow separation from the glottal wall which provides sufficient negative damping to maintain oscillatory motion of the folds. The separation zone changes during the inward and outward oscillations of the folds in the manner illustrated in Figure 4, which shows a cross-section of the flow in a median plane normal to the flow and parallel to the short sides of the rectangular model glottis of Figure 3b, with the flow through the glottis being from left to right. In condition (a) dζ/dt > 0 and the glottis is expanding: the incoming flow separates at the local ‘trailing edge’ at the front edge A of the glottis, so that the pressure applied over the glottal walls is uniform and approximately the same as the uniform free space pressure. In case (b) dζ/dt ≤ 0, the glottal area is decreasing and vorticity is now shed from the neighbourhood of the trailing edge at B. In voiced speech the effective Strouhal number of the unsteady motion within the glottis fo[ell]g/V ~ 0.1 (where the mean flow speed across the glottis of axial extent [ell]g is typically of order V ~ 20 m/s), which is small enough for it to be assumed that, over a cycle of oscillation at frequency fo, the fluid passes through a continuous sequence of quasi-static states of motion.

Figure 3
(a) Large aspect ratio elliptic approximation to the shape of the glottal cross-section during phonation; (b) Rectangular approximation used in the text.
Figure 4
Operating stages of the ‘one mass’ model of the glottis:

During expansion vorticity is shed from the fore edge A of the glottis in the manner illustrated in Figure 4a, resulting in the formation of a jet with a contraction ratio (relative to the glottal cross-section) σ ~ 0.62; when the glottis begins to contract the point of free shear layer formation rapidly jumps to the aft edge B of the glottis (Figure 4b), beyond which the flow expands into a jet whose cross-section is marginally larger than that of the orifice, with σ ~ 1.15. The resulting attached accelerated flow within the glottis produces a large decrease in the glottal wall pressure (below that in the downstream region) and correspondingly large ‘suction’ forces that reinforce the rate of decrease in glottal area. Further details of the calculation of these forces are given by Howe and McGowan (2009).

Symmetric oscillatory fluctuations in the glottal width 2ζ(t) of the model of Figure 4 are assumed to be governed by the following oscillator equation:

d2ζdt2+γodζdt+Ω2(ζζ0)=F(t)m,
(2.1)

where m is the effective mass of each fold, ζ0 is the equilibrium semi-width of the glottis, γo is a structural damping coefficient, and Ω is the radian resonance frequency determined by the elastic properties of the fold. The aerodynamic force F is the suction force exerted on the glottal wall when dζ/dt ≤ 0 (Case (b) of Figure 4) and is otherwise zero, i.e. (see Howe and McGowan (2009) for further details)

F={1332πρoUσ2σΔ,dζdt0,0,dζdt>0,
(2.2)

where Δ [equivalent] Δ(t) = 4(t) is the cross-sectional area of the glottis and Uσ [equivalent] Uσ (t) is the contracted jet velocity in the uniform downstream section of the jet (cf. Figure 5).

Figure 5
Streamsurfaces of the potential function Φ(y, τ) intersecting the jet boundary with outward normal n.

The suction force formula (2.2) is based on a precise modelling of the hydrodynamic flow within and near the glottis using ‘free streamline theory’ (unlike, for example, the heuristic negative damping mechanism proposed for the maintenance of self-sustained oscillations by Fulcher et al. (2006)). Its discontinuous behaviour arises because of the neglect in the mathematical model for F of the finite streamwise length [ell]g of the glottis. This implies that the time ~ [ell]g/V required for the separation point to move between the ends of the glottis is negligible, so that the adjustment occurs instantaneously. Provided this limitation is understood the reader should experience no difficulty in interpreting any consequential (and small) discontinuities in predicted acoustic waveforms. Such discontinuities were avoided by Zanartu et al. (2007) by using an arbitrary smoothing of discontinuities associated with abrupt changes in the ‘glottis discharge coefficient’.

Similarly, according to this model mutual impact of the vocal folds occurs when ζ → 0. In solving equation (2.1) it is assumed that this impact is inelastic, so that dζ/dt is then reduced to zero.

3. The aeroacoustics of voicing

3.1 The steady lung pressure

The ultimate source of sound is the uniform volume flow from the planar volume source spanning the lung plenum at x = −[l with macron] (Figure 1). The idealized subglottal system has acoustically rigid walls and it is assumed that all relevant frequencies, including the resonance frequencies of the subglottal system, are small enough that only plane waves can propagate. Similarly, dissipation within the flow is ignored, except in so far as acoustic energy is lost by radiation into free space through the glottis, and by radiation into the lung plenum from the subglottal tube. The motion produced by the source in a homogeneous fluid may therefore be regarded as irrotational and defined by a velocity potential ϕ. If the flow starts at time t = τ0, say, the potential is determined by

(1co22t22)ϕ=qH(tτ0)δ(x+¯),
(3.1)

where co is the speed of sound. The volume source on the right hand side is assumed to have the constant strength q per unit area in the plane x = −[l with macron] and H(tτo) is the Heaviside step function, equal to 0, 1 according as t [less, greater] τ0. It is convenient to choose the source time τ0 = −[l with macron]/co. Then initially the solution of equation (3.1) is a plane wave ϕ0 radiating from the source towards the subglottal tract, given by

ϕ0=coq2(txco)H(txco),¯<x<.
(3.2)

This wave sets up a steady flow at speed [partial differential]ϕ0/[partial differential]x = q/2 from the source towards the subglottal tube.

A plane wave ϕR is reflected back into the lung plenum when this wave reaches the area change at the junction at x = −[ell]. It is given by

ϕR=coq2(ALAAL+A)(t+xco+2co)H(t+xco+2co),x<,
(3.3)

where (An external file that holds a picture, illustration, etc.
Object name is nihms148284ig2.jpgAn external file that holds a picture, illustration, etc.
Object name is nihms148284ig1.jpg)/(An external file that holds a picture, illustration, etc.
Object name is nihms148284ig2.jpg + An external file that holds a picture, illustration, etc.
Object name is nihms148284ig1.jpg) is the junction reflection coefficient (Pierce 1989). The resulting initial pressure established just to the left of the junction is therefore

p=ρot(ϕ0+ϕR)=ρocoqALAL+A
(3.4)

where ρo is the mean fluid density.

We now define the pressure pI by

pI=ρocoqALAL+A,sothatq=pIρoco(AL+AAL).
(3.5)

The pressure pI is the steady pressure produced by the volume source at the remote end x = −[ell] of the subglottal tube in the absence of sound generation and reflections at the glottis. By continuity, the pressure wave transmitted into the subglottal tube when ϕ0 arrives at the junction also has magnitude pI, and sound begins to be generated by the glottis at x ~ 0 after the arrival of the pressure wave

p=pIH(txco)
(3.6)

at time t = 0.

In the following we shall assume that the pressure pI is prescribed (typically of the order of one kPa in magnitude).

3.2 Formal representation of sound radiated from the glottis

In the absence of dissipation entropy changes within the flow can be neglected and sound production is governed by a simplified form of the equation of vortex sound (Howe 1998). This determines the acoustic field in terms of the fluctuations in the total enthalpy B, given in homentropic flow by

B=dpρ+12v2,
(3.7)

where ρ= ρ(p) is the local density and v is the velocity. In those regions where the flow is irrotational Bernoulli’s equation implies that B [equivalent][partial differential]ϕ/[partial differential]t. It then follows by combination of the vortex sound equation (Howe 1998) and the time differential of equation (3.1), that the composite equation governing sound production by the source flow q and the jet flow interaction with the glottis can be cast in the form

(1co22t22)B=pIρoco(AL+AAL)δ(t+¯co)δ(x+¯)+div(ωv).
(3.8)

The first term on the right is the contribution from the flow source q of equation (3.1), with q replaced by its definition (3.5) in terms of pI and the start time τ0 of equation (3.1) replaced by −[l with macron]/co. The divergence source term involves the vorticity ω = curl v and is confined principally to the free shear layers of the jet.

Equation (3.8) is solved by use of Green’s function G(x, y, t, τ), which corresponds to the sound produced by an impulsive point source at y = (x, y, z′) at time τ, and satisfies

(1co22t22)G=δ(xy)δ(tτ),G=0fort<τ.
(3.9)

The normal derivative [partial differential]G/[partial differential]xn = 0 vanishes on the rigid walls of the vocal tract, and is also required to vanish on the instantaneous position of the moving walls of the glottis.

The application of Green’s theorem and the radiation condition in the usual way (Baker and Copson 1969; Landau and Lifshitz 1987; Crighton et al. 1992; Howe 1998) then supplies the following formal representation of the solution of equation (3.8)

B(x,t)=pIρoco(AL+AAL)SLG(x,(¯,y,z),t,¯/co)dydz(ωv)(y,τ)·Gy(x,y,t,τ)d3ydτ+S(τ)G(x,y,t,τ)vτ(y,τ)·dS(y)dτ,
(3.10)

where in the linear acoustic region B(x, t) [equivalent][partial differential]ϕ/[partial differential]tp(x, t)o.

The first term on the right hand side represents the direct sound generated by the lung source q at x = −[l with macron], the integration being over the the cross-section SL of the lung plenum occupied by this source. The second integral gives the sound generated by vortex sources in the unsteady jet emerging from the glottis, the integration being over the region of space wherein ω ≠ 0 and over all source times −∞ < τ < ∞. The final term involves a surface integration of the normal component of velocity over the walls S(τ) of the glottis, and accounts for the production of ‘monopole’ sound by the volumetric fluctuations of the oscillating vocal folds. This mechanism is known to make a relatively small contribution to voicing (see, e.g., Zhao et al. 2002; Howe and McGowan 2007) and will henceforth be ignored.

3.3 Green’s function

There is no closed form representation of Green’s function that is applicable under all conditions. But for the present purposes it is only necessary to be able to calculate the sound using equation (3.10) in the limit in which the dominant acoustic wavelengths are large compared to the diameter of the vocal tract. In these circumstances only plane waves can propagate within the subglottal tube and in the lung plenum x < −[l with macron], and approximate representations for G(x, y, t, τ) are easily derived. The details of the calculations are summarised in the Appendix.

To fix ideas we consider the calculation of the sound in free space at the point x of Figure 1. The two principal source terms in the integrals of equation (3.10) are situated (i) at points y on the volume source plane x = −[l with macron] within the lungs, i.e. within the region labelled L in Figure A.1 of the appendix, and (ii) at points y in the immediate neighbourhood of the glottis (region E of Figure A.1). It is shown in the appendix that the corresponding approximate functional forms of Green’s function are GL and GE given by:

Figure A1
The regions L, S, E, A used to define the approximate Green’s function.

Case (i), x in free space, y [equivalent] (x, y, z′) in the lung plenum x′ < −[ell]:

GL(x,y,t,τ)=T2πxn=0Rnδ(tτ+xco[x+2nL+E]co);
(3.11a)

Case (ii), x in free space, y [equivalent] (x, y, z′) in the neighbourhood of the glottis:

GE(x,y,t,τ)=12πxδ(tτxco)+Φ(y,τ)2πcox{n=0Rnδ(tτ[x+2nL+E]co)+n=0Rn+1δ(tτ[x+2(nL+)+E]co)}.
(3.11b)

In these formulae the prime on the delta functions denotes differentiation with respect to the argument and

R=(ALAAL+A),T=2AAL+A
(3.12)

are respectively the long wavelength reflection and transmission coefficients for a plane wave incident on the junction at x = −[ell] from within the subglottal tube. The function Φ(y, τ) is a velocity potential describing flow through the glottis in its state at time τ; it is a solution of Laplace’s equation with vanishing normal derivative on all boundaries, and it is normalized such that

Φ(y,τ)A2πyasyinfreespace,xatlargedistancesfromtheglottisinthesubglottaltract,}
(3.13)

where (τ)O(A) is the glottis ‘end correction’ (Rayleigh 1945; Howe 1998). It is shown in the Appendix that in calculating the phase of successive terms in the expansions (3.11a, 3.11b) it is permissible to replace all occurrences of [ell]′ by a mean value

E=,
(3.14)

where the angle brackets denote a suitable average over a cycle of operation of the glottis. The length L is just the length of the subglottal tube augmented by the end correction: L = [ell] + [ell]E.

The end correction [ell]′ becomes very large when the cross-section of the glottis becomes small. This can cause the effective acoustic length of the subglottal tract to vary during a period of phonation, producing a variation of subglottal resonance frequencies. However near closure, when [ell]′ is large, the acoustic amplitude or volume flux through the glottis is small. Thus, for those times when the acoustic flux is significant we can use a smaller, typical value of [ell]E. Numerical results given by Howe and McGowan (2007), for example, suggest that, at least for the purposes of the present calculation, [ell]E/h ≤ 5.

3.4 The free space radiation

We are now ready to evaluate the sound radiated from the glottis, given by equation (3.10). For the first integral on the right hand side Green’s function must be taken in the form GL of equation (3.11a). Taking account of the definition (3.12) of An external file that holds a picture, illustration, etc.
Object name is nihms148284ig3.jpg we find

pIρoco(AL+AAL)SLG(x,¯,y,z,t,¯/co)dydz=A2πxt{2pIρocon=0RnH(t[x+2nL+E]co)}.
(3.15)

For the second, vortex integral on the right of (3.10) we use the form G = GE, for which it is necessary to evaluate

(ωv)(y,τ)·Φy(x,τ)d3y.

To do this we use the method of Howe and McGowan (2007) and observe that in a first approximation the vorticity is confined to the free streamsurfaces at the edges of the jet emerging from the glottis, as illustrated in Figure 5. This figure also displays the instantaneous hypothetical family of streamsurfaces of the flow through the glottis determined by the potential Φ (y, τ). On the jet interface vorticity is convected at half the local jet velocity, and we can put ωv=12Uσ2δ(s)n, where s[perpendicular] is distance measured in the direction of the outward unit normal n of the jet and Uσ is the jet velocity just inside the shear layer. The main contribution to the integral is from that section of the jet very close to the glottis where the jet interface is cut by the Φ-streamlines. The velocity Uσ does not vary significantly over this short distance along the jet, and we therefore have

(ωv)(y,τ)·Φy(x,τ)d3yUσ2(τ)2SJΦy(y,τ)·dSAUσ22,
(3.16)

where the second integral is over the curved surface SJ of the jet, and we have used the normalisation condition (3.13) of Φ (y, τ).

Hence

(ωv)(y,τ)·Gy(x,y,t,τ)d3ydτ=A2πcoxt{12Uσ2(t(x+E)co)+n=112Rn[Uσ2(t(x+2nL+E)co)+Uσ2(t(x+2nLE)co)]}
(3.17)

When the results (3.15) and (3.17) are substituted into the right hand side of (3.10) (with the final monopole integral discarded), and B(x, t) is replaced by p(x, t)o in the acoustic domain, we find that the free space far field sound is given by

p(x,t)ρoA2πcoxt{2pIρon=0RnH(t[x+2nL+E]co)12Uσ2(t(x+E)co)n=112Rn[Uσ2(t(x+2nL+E)co)+Uσ2(t(x+2nLE)co)]},x.
(3.18)

4. The glottal volume flux

4.1 The volume flux equation

Equation (3.18) defines an omnidirectional acoustic field radiating into the half-space x > 0 with all of the characteristics of a monopole. As such the sound can be expressed in terms of the unsteady volume flux Q = ΔσUσ, where Δσ is the asymptotic contracted jet cross-sectional area downstream of the glottis.

Let ϕ(x, t) denote the velocity potential in the acoustic far field, where p = −ρo[partial differential]ϕ/[partial differential]t. Then equation (3.18) is equivalent to

ϕ(x,t)A2πcox{2pIρon=0RnH(t[x+2nL+E]co)12Uσ2(t(x+E)co)n=112Rn[Uσ2(t(x+2nL+E)co)+Uσ2(t(x+2nLE)co)]},x.
(4.1)

But, we can also write

ϕ(x,t)Q2πx(txco)Δσ2πxUσ(txco),x.
(4.2)

Therefore, if we now define

Mσ(t)=Uσ(t)co,
(4.3)

then equations (4.1), (4.2) imply that

Mσ2(t)+2ΔσAMσ(t)=4pIρoco2n=0RnH(t2nLco)n=1Rn{Mσ2(t2nLco)+Mσ2(t2nLco+2Eco)}.
(4.4)

This equation determines the glottal volume flux

Q(t)=ΔσcoMσ(t),
(4.5)

provided Δ [equivalent] Δ(t) is also known. Equation (4.4) must therefore be solved simultaneously with equation (2.1) for ζ= Δ/4a.

To do this it is convenient in the general case to express Mσ (t) directly in terms of pI and the values of Mσ at earlier times, as follows

Mσ(t)={(σΔ(t)A)2+4pIρoco2n=0RnH(t2nLco)n=1Rn[Mσ2(t2nLco)+Mσ2(t2nLco+2Eco)]}12σΔ(t)A.
(4.6)

Here we have taken the larger (and generally positive) root of the quadratic equation (4.4). This is in accord with experiment that requires the net volume flow to be out of the glottis, and also represents the continuous solution of the problem starting from t = 0.

At large times, when the voiced radiation assumes an invariant waveform we can use the asymptotic formula

n=0RnH(t2nLco)=AL+A2A,t.

If we now put

PI=(AL+A2A)pI,
(4.7)

then (4.6) becomes

Mσ(t)=(σΔ(t)A)2+4PIρoco2n=1Rn[Mσ2(t2nLco)+Mσ2(t2nLco+2Eco)]σΔ(t)A,t.
(4.8)

We can check the consistency of this result in the irrotational limit, in the absence of jetting from the glottis. Nonlinear terms are then negligible and the glottal outflow expands irrotationally into the half-space x > 0, i.e. σ → ∞. Then Mσ2 can be neglected in the square root on the right of (4.8) which, because (σΔ/A)24PI/ρoco2, can be expanded to first order in PI yielding (with the aid of equations (4.7) and (3.5))

Q=σcoΔ(t)Mσ(t)=2APIρoco=(AL+A)ρocopIALq.
(4.9)

This confirms that, following the passage of all transients and when glottal nonlinearity is unimportant, all of the volume flow An external file that holds a picture, illustration, etc.
Object name is nihms148284ig2.jpgq from the lung plenum volume source flows out steadily through the glottis.

4.2 Resonant oscillations

When [ell]E [double less-than sign] L the glottal oscillations occur without impeding significantly the volume flow from the lungs. The resonance oscillations then correspond to standing waves in the subglottal tract in which the period 1/f0 of the motion of the vocal folds is just equal to the travel time 2L/co required by sound waves for one ‘round trip’ between the opposite ends of the tract. In that case, and for large times, we have

Mσ(t2nLco)=Mσ(t),Mσ(t2nLco+2Eco)=Mσ(t2co).

Equation (4.4) then becomes

Mσ2(t)+2ΔσTAMσ(t)=4TPIρoco2RMσ2(t2co),
(4.10)

so that

Mσ(t)=(σΔ(t)TA)2+4TPIρoco2RMσ2(t2co)σΔ(t)TA,t.
(4.11)

Similarly, when fo[ell]E/co [double less-than sign] 1 we have the corresponding further simplifications

Mσ2(t)+2ΔσALMσ(t)=4PIρoco2AAL,
(4.12)

and

Mσ(t)=AAL{(σΔ(t)A)2+4PIρoco2ALAσΔ(t)A},t.
(4.13)

5. Numerical results

5.1 Parameter values

Previous numerical and analytical modelling (e.g. Zhao et al. 2002; Zanartu et al. 2007; Howe and McGowan 2007, 2009) indicates that the predicted waveform for steady excitation of the glottis by a nominally constant over pressure pI arriving from the lungs rapidly attains a steady periodic form. In the present case predictions of the unsteady volume flux Q(t) = ΔσcoMσ (t) through the glottis are governed by the simultaneous solution of the vocal fold equation (2.1) for ζ(t) and the appropriate equation from Section 4 that determines Mσ(t) in terms of pI and Δ [equivalent] 4(t).

The numerical results discussed in this Section illustrate predictions when the natural frequency of the folds f0 [equivalent] Ω/2π = 125 Hz, typical of an adult male, and for pI = 10 cm of water (~ 1 kPa) when co = 340 m/s, ρo = 1.23 kg/m3. The following values have also been assumed for other parameters entering the formulae:

h=10mm,3/h=1.3,m=0.5×104kg,γo/f0=0.1.
(5.1)

We have used values of the jet contraction ratio recommended by Howe and McGowan (2009)

σ={0.62,fordζ/dt>0,1.15,fordζ/dt<0.
(5.2)

The folds are first set into motion from a state of rest by the arrival at the glottis of the pressure pI, because the initial flow through the glottis is irrotational and it necessarily produces a non-zero value for the suction force F (t) on the right of equation (2.1). The initial rest position of the folds is taken to be ζ= ζ0, with dζ/dt = 0, where ζ0/h = 0.05, and the solution is advanced by means of a fourth order Runge-Kutta procedure.

5.2 The infinitely long subglottal tube

The dissipation of acoustic waves radiated towards the lungs may be assumed to suppress reflections from the inner end x = −[ell] of the subglottal tract when [ell] becomes very large. This is mathematically equivalent to setting R = 0 in the formulae of Section 4, in which case equation (4.6) reduces to

Mσ(t)=(σΔ(t)A)2+4pIρ0co2σΔ(t)A,t.
(5.3)

Figure 6a displays two cycles of the corresponding volume flux waveform

Figure 6
(a) Two cycles of the volume flux through the glottis in the limiting case of an infinitely long subglottal tract (R = 0) when f0 = 125 Hz, ζ0/h = 0.05 and pI = 10 cm of water.

QQoρoco2pIσ(t)Mσ(t)ζ(t)h
(5.4)

plotted as a function of the nondimensional time f0t, where Qo = An external file that holds a picture, illustration, etc.
Object name is nihms148284ig1.jpgpIoco is the volume flux of the initial pressure pulse incident on the glottis. Also shown in Figure 6b is the calculated variation of the semi-width ζ(t)/h of the glottis.

The rapid variation in Q near f0t=12,32 is caused by the abrupt increase in the value of the jet contraction ratio σ at those times at which dζ/dt becomes negative. The near discontinuity occurs because of the correspondingly rapid motion of the jet separation point (Figure 4) over the small time interval ~ [ell]g/Uσ, which is neglected in the present calculation. This is evident from Figure 7, which shows predicted uniformly smooth variations of Q/Qo obtaining when the same calculation is performed for a constant mean value of the contraction ratio σ = 0.8.

Figure 7
(a) Two cycles of the volume flux through the glottis in the limiting case of an infinitely long subglottal tract (R = 0) and for a fixed value σ = 0.8 of the jet contraction ratio, when f0 = 125 Hz, ζ0/h = 0.05 and pI = 10 cm ...

5.3 Influence of the length of subglottal tube for lightly damped waves

Two damping mechanisms are included in the present theory: the internal damping of the vocal fold material, specified by the coefficient γo of equation (2.1), and damping due to radiation into the lung plenum, determined by the reflection coefficient R. Light damping by the lungs accordingly corresponds to |R| ~ 1; it is then likely that the multiple reflection of sound from opposite ends of the subglottal tube will have a significant impact on the volume flux through the glottis. The magnitude of this effect must also depend on the length [ell] of the subglottal tube, or rather on the effective acoustic length L = [ell] + [ell]E.

Figure 8 depicts the calculated variation of Q/Qo for one cycle when R = 0.9 for a range of increasing subglottal lengths L in the limit in which the end correction [ell]E of the glottis is neglected (or, more strictly, when f0[ell]E/co is sufficiently small that the approximation [ell]E = 0 is assumed to be valid in equation (4.6), wherein the series expansions were truncated at n = 50). Case (a) of the figure, where 2f0L/co = 0.2, would correspond roughly to the situation in an adult male if we had taken [ell]E ~ 10 cm (with [ell] ≈ 17 cm). Figures 8b–8e show the progressive changes in the volume flux waveform predicted to occur in an experimental simulation in which the length of the subglottal tube is increased in stages to a maximum at which 2f0L/co = 1. In case (e) the duct length is one half of the acoustic wavelength at frequency f0 (i.e. f0 then corresponds to the resonance frequency of an open ended tube of length L) – the waveform is seen to take its maximum value Q/Qo ~ 2 in this resonant case, and to resemble in shape the profile for an infinite duct (R = 0, Figure 6). Peak values are smaller for the other cases shown in the figure. However, Cases (b), (c) indicate a strong influence of the second harmonic, occurring when the duct length corresponds approximately to a quarter-wave resonator.

Figure 8
Dependence of the volume flux waveform on the length L of the subglottal duct when |R| = 0.9, f0 = 125 Hz, 2f0[ell]E/co ~ 0, ζ0/h = 0.05 and pI = 10 cm of water; for 2f0L/c0 = (a) 0.2, (b) 0.45, (c) 0.6, (d) 0.8, (e) 1.0.

The calculations have been repeated (see Figure 9) for the same conditions, but the retention of [ell]E in equation (4.6) for the case [ell]E = 13 cm (so that L = [ell] + [ell]E ~ 30 cm). Increasing the mean end-correction of the glottis is equivalent to increasing its acoustic resistance. In the notation of the appendix, this value of [ell]E corresponds to ko[ell]E ~ 0.3, or 2f0[ell]E/co = 0.1, (where ko [equivalent] 2πf0/co is the acoustic wavenumber) which is the largest value for which (ko[ell]E)2 can be said to be small compared to unity, a condition required for the validity of the Green’s function used in these calculations.

Figure 9
Dependence of the volume flux waveform on the length L of the subglottal duct when |R| = 0.9, f0 = 125 Hz, 2f0[ell]E/co = 0.1, ζ0/h = 0.05 and pI = 10 cm of water; for 2f0L/c0 = (a) 0.2, (b) 0.45, (c) 0.6, (d) 0.8, (e) 1.0.

A comparison of corresponding plots in Figure 8 (2f0[ell]E/co ≈ 0) and Figure 9 (2f0[ell]E/co = 0.1, [ell]E = 13 cm) indicates very similar volume flux profiles and amplitudes, with the second harmonic appearing in both cases when 2f0L/co ~ 0.5.

5.4 Heavily damped waves

The possible influences of subglottal wave motions on the glottal flow are much reduced when the damping increases, which is achieved theoretically by reducing the magnitude of the reflection coefficient R. Volume flux profiles are plotted in Figure 10 for cases where |R| = 0.5 and [ell]E = 13 cm (2f0[ell]E/co = 0.1). Evidently, there are now only relatively minor changes in the waveform as 2foL/co increases from 0.2 in Case (a) to the resonance configuration of Case (e), with a very marginal increase in amplitude at resonance.

Figure 10
Dependence of the volume flux waveform on the length L of the subglottal duct when |R| = 0.5, f0 = 125 Hz, 2f0[ell]E/co = 0.1, ζ0/h = 0.05 and pI = 10 cm of water; for 2f0L/c0 = (a) 0.2, (b) 0.45, (c) 0.6, (d) 0.8, (e) 1.0.

6. Conclusion

Experimental laboratory models designed to replicate the interaction between flow in the subglottal tube and the glottis can be influenced by sound waves in the subglottal system that are excited by the interaction. It is therefore important to take proper account of the acoustic properties of the whole system in order to obtain a proper interpretation of experimental results. The theory of this paper is roughly applicable to those studies of the subglottal region where unwanted acoustic feedback from the supraglottal vocal tract needs to be suppressed, and which is achieved experimentally by allowing the glottis to communicate directly with the ambient free space.

Our nonlinear mathematical model of these interactions is applicable in phonation instances where the glottal cross-section is ‘large’ over most of a typical cycle, so that its acoustical ‘end correction’ does not become too large. Numerical predictions indicate that a change in the length of the experimental subglottal tube can produce alterations in the unsteady (sound generating) glottal volume flux during excitation by a nominally steady volume flow of air from the lungs. These modifications become very significant in a lightly damped experimental system, where proper account is not taken of absorption damping of the sound by the lung complex. There is a noticeable but small increase in volume flux magnitude at the resonance occurring when the effective acoustic length L of the subglottal tube is equal to half the acoustic wavelength; perhaps of more significance is the strong appearance of second harmonic generation at intermediate lengths, when L is roughly one quarter of the acoustic wavelength. On the other hand, changes in the volume flux with the length L tend to be relatively small when lung damping is large.

These detailed conclusions must necessarily be dependent to some extent on the mechanical model of the glottis. The discontinuous nature of the aerodynamic edge force on the glottis is of no importance – this is a common feature of other models that is usually suppressed without mention. The prediction that nonlinear coupling between the glottal motion and the subglottal fluid and acoustic motions produces a modified volume flux through the glottis appears to be quite general.

Acknowledgments

This work was supported by a subaward of grant No. 1R01 DC009229 from the National Institute on Deafness and other Communication Disorders to the University of California, Los Angeles.

Appendix: Approximate Green’s function

The solution of equation (3.9) is required at points x in the free space in region A of Figure A1 at distances |x| [dbl greater-than sign] 2a from the glottis when y is in region E close to the glottis.

Put

G(x,y,t,τ)=12πG¯(x,y;ω)eiω(tτ)dω,
(A.1)

so that G satisfies the inhomogeneous Helmholtz equation (2+ko2)G¯=δ(xy), ko = ω/co. The characteristic frequency is taken to be sufficiently small ( koA1) that only plane waves can propagate within the subglottal and lung ducts.

The functional form of G(x, y; Ω) is derived using reciprocal relation G(x, y; ω) = G(y, x; ω) (Rayleigh 1945), by placing the point source at x in free space and solving for G(y, x; ω) as a function of y in the vicinity of the glottis (cf. Howe 2003, 2005). The source generates a spherical wave incident on the glottis; the spherical wave is reflected from the glottis flange, so that the net exterior field at E in Figure A1 can be taken in the form

G¯eikox2πx+γeikoy2πy,
(A.2)

where the second term on the right (involving an amplitude γ to be determined) represents a spherical wave radiated from the glottis in response to the forcing of flow through the glottis by the incident wave and by waves reflected back and forth within the subglottal tube.

Within and in the immediate vicinity of the glottis we write

Gα+βΦ(y),
(A.3)

where α, β are constants to be determined, and Φ*(y) is defined with the asymptotic properties (3.13) for some particular configuration of the glottis, which for the moment is assumed to be of fixed shape. This expression for G must merge into a combination of plane waves within the body S of the subglottal tract, where

G¯ASeikox+BSeikox,<x<0,
(A.4)

for suitable constants AS, BS. In the lung region L of Figure A1, waves entering from the subglottal tube are absorbed (nominally at x = −∞) in the lung, where we therefore put

G¯CLeikox,x<,
(A.5)

for some constant CL.

The coefficients α, β,γ, AS, BS, CL in (A.2) – (A.5) are determined by matching the different representations of G where neighbouring regions of validity overlap. At the junction between the lung and subglottal tube it is sufficient to equate the potentials and volume flux on either side of x′ = −[ell], which yields

AS=(ALA2A)CLe2iko,BS=(AL+A2A)CL
(A.6)

Similarly equating expressions for G and [partial differential]G/[partial differential]x′ just to the left of the glottis in Figure A1 (where ko|x′| [double less-than sign] 1 and the second approximation of (3.13) is applicable), we find

αβ={(AL+A)(ALA)eiko}CL2A,β=iko{(AL+A)+(ALA)eiko}CL2A.}
(A.7)

Finally, outside the glottis in the overlap of regions E and A, the first of conditions (3.13) yields

α=eikox2πx+ikoγ2π,β=γA.
(A.8)

Equations (A.6)(A.8) must now be solved for the coefficients in terms of the incident spherical wave from x. However, because the sources in equation (3.8) are either within the lung or in the jet just outside the glottis, only the values of α, β and CL are actually needed. In particular, noting that ko[ell]′ and ko2A are both small, we find

β=iko(1+Re2iko)eiko(x+)ko2A/2π2πx(1Re2ikoLko2A/π),
(A.9)

where L = [ell] + [ell]′ and R is defined as in equation (3.12).

The factor eko2A/2π in this formula represents the action of dissipation produced by radiation into free space (via the second spherical wave in (A.2)). Because ko2A1 this dissipation is small compared to the absorption by the lungs, whose magnitude is determined by that of |R| < 1, and it will therefore be neglected, leading to

β=iko(1+Re2iko)eiko(x+)2πx(1Re2ikoL).
(A.10)

To the same order of approximation it is found that

α=eikox2πx,CL=Teiko(x+)2πx(1Re2ikoL),
(A.11)

where An external file that holds a picture, illustration, etc.
Object name is nihms148284ig3.jpg is defined in (3.12).

The space-time form of Green’s function for source points y in regions E and L are now obtained by substitution of the formulae (A.3) and (A.5) into the integral (A.1). To illustrate the procedure consider the evaluation of

GL(x,y,t,τ)=12πCLeiω(tτ+x/co)dω.
(A.12)

The condition |R| < 1 permits CL to be expanded in the form

CL=T2πxn=0Rne2inkoL+iko(x+).
(A.13)

Substitution into (A.12) and term by term integration then yields the result (3.11a) of the main text.

A similar expansion is used to obtain the formal representation (3.11b) for GE. Note, however, that the above derivation by Fourier transformation is strictly applicable for a glottis of fixed configuration, defined by time-independent values of Φ*(y) and [ell]′. But each of the expanded representations (3.11a, b) clearly identifies an infinite sequence of separate interactions with the glottis occurring at distinct source times. The formulae can therefore be generalized to accommodate slow variations in glottal geometry simply by associating each interaction at the appropriate source time τ with the corresponding potential function, such that Φ* → Φ (y, τ) and with the values of the [ell]′s chosen with the appropriate time delay determined by the number of internal reflections in the subglottal tract prior to emission through the glottis and radiation to the far field region A of Figure A.1. This is the procedure successfully applied and described in more detail by Howe and McGowan (2007).

However, this means that there is a succession of wave elements at x in the far field, corresponding to the successive terms in the expansion (3.11a), each of which has a complicated phase dependence. Thus, the phase of the nth component of (3.11a) would depend on n different values of [ell]′ determined by the retarded times at which internal reflection occurs at the glottis prior to radiation into free space. When |R| is large many terms in the series expansion (3.11a) are required, and many multiple reflections within the tract are important. In this case the relevant values of [ell]′ will tend to be distributed over a large number of realizations of the state of the glottis during a typical cycle, and to a good approximation the phase can be calculated by replacing each occurrence of the end-correction [ell]′ by its mean value [ell]E defined as in (3.14). The same can be done when |R| is small, because in this case the number of relevant reflections is small and the corresponding overall phase shift introduced by the reflections is small.

Footnotes

Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Contributor Information

M. S. Howe, Boston University, College of Engineering, 110 Cummington Street, Boston MA 02215.

R. S. McGowan, CReSS LLC, 1 Seaborn Place, Lexington MA 02420.

References

  • Alipour F, Scherer RC. Effects of oscillation of a mechanical hemilarynx model on mean transglottal pressures and flows. Journal of the Acoustical Society of America. 2001;110:1562–1569. [PubMed]
  • Austin SF, Titze IR. The effect of subglottal resonance upon the vocal fold vibration. Journal of Voice. 1997;11:391–402. [PubMed]
  • Baker BB, Copson ET. The Mathematical Theory of Huygens’ Principle. 2. Oxford University Press; 1969.
  • Crighton DG, Dowling AP, Ffowcs Williams JE, Heckl M, Leppington FG. Modern Methods in Analytical Acoustics (Lecture Notes) Springer-Verlag; London: 1992.
  • Fant G. Acoustic Theory of Speech Production. Mouton: The Hague; 1960.
  • Flanagan JL. Some properties of the glottal sound source. Journal of Speech, and Hearing Research. 1958;9:99–116. [PubMed]
  • Flanagan JL, Landgraf LL. Self-oscillating sources for vocal-tract synthesizer. IEEE Transactions Audio and Electroacoustics. 1968;AU-16:57–64.
  • Flanagan JL. Speech Analysis Synthesis and Perception. 2. Springer-Verlag; New York: 1972.
  • Fulcher LP, Scherer RC, Melnykov A, Gateva V, Limes ME. Negative Coulomb damping, limit cycles, and self-oscillation of the vocal folds. American Journal of Physics. 2006;74:386–393.
  • Gupta V, Wilson TA, Beavers GS. A model for vocal cord excitation. Journal of the Acoustical Society of America. 1973;54:1607–1617. [PubMed]
  • Howe MS. Acoustics of Fluid-Structure Interactions. Cambridge University Press; 1998.
  • Howe MS. On the infrasound generated when a train enters a tunnel. Journal of Fluids and Structures. 2003;17:629–642.
  • Howe MS. The genetically optimized tunnel-entrance hood. Journal of Fluids and Structures. 2005;23:1231–1250.
  • Howe MS, McGowan RS. Sound generated by aerodynamic sources near a deformable body, with application to voiced speech. Journal of Fluid Mechanics. 2007;592:367–392.
  • Howe MS, McGowan RS. On the single-mass model of the vocal folds. Fluid Dynamics Research. 2009 in press. [PMC free article] [PubMed]
  • Ishizaka K, Matsudaira M, Kaneko T. Input acoustic-impedance measurement of the subglottal system. Journal of the Acoustical Society of America. 1976;60:190–196. [PubMed]
  • Joliveau E, Smith J, Wolfe J. Vocal tract resonances in singing: The soprano voice. Journal of the Acoustical Society of America. 2004;116:2434–2439. [PubMed]
  • Landau LD, Lifshitz EM. Fluid Mechanics. 2. Pergamon; Oxford: 1987.
  • Pierce AD. Acoustics, An introduction to its principles and applications. American Institute of Physics; 1989.
  • Rayleigh Lord. Theory of Sound. Vol. 2. Dover; New York: 1945.
  • Thomson SL, Mongeau L, Frankel SH. Aerodynamic transfer of energy to the vocal folds. Journal of the Acoustical Society of America. 2005;118:1689–1700. [PubMed]
  • Titze IR. The physics of small-amplitude oscillation of the vocal folds. Journal of the Acoustical Society of America. 1988;83:1536–1552. [PubMed]
  • Titze IR. Nonlinear source-filter coupling in phonation: Theory. Journal of the Acoustical Society of America. 2008;123:2733–2749. [PubMed]
  • Titze IR, Story BH. Acoustic interactions of the voice source with the lower vocal tract. Journal of the Acoustical Society of America. 1995;101:2234–2243. [PubMed]
  • Zanartu M, Mongeau L, Wodicka GR. Influence of acoustic loading on an effective single mass model of the vocal folds. Journal of the Acoustical Society of America. 2007;121:1119–1129. [PubMed]
  • Zhang Z, Neubauer J, Berry DA. The influence of subglottal acoustics on laboratory models of phonation. Journal of the Acoustical Society of America. 2006a;120:1558–1569. [PubMed]
  • Zhang Z, Neubauer J, Berry DA. Aerodynamically and acoustically driven modes of vibration in a physical model of the vocal folds. Journal of the Acoustical Society of America. 2006b;120:2841–2849. [PubMed]
  • Zhang Z, Neubauer J, Berry DA. Influence of vocal fold stiffness and acoustic loading on flow-induced vibration of a single-layer vocal fold model. Journal of Sound and Vibration. 2009;322:299–313. [PMC free article] [PubMed]
  • Zhao W, Zhang C, Frankel SH, Mongeau L. Computational aeroacoustics of phonation, Part I: Computational methods and sound generation mechanisms. Journal of the Acoustical Society of America. 2002;112:2134–2146. [PubMed]