Direct-numerical simulation of the glottal jet and vocal-fold dynamics in a three-dimensional laryngeal model

Direct-numerical simulation of the glottal jet and vocal-fold dynamics in a three-dimensional laryngeal model X. Zheng, R. Mittala) and Q. Xue Department of Mechanical Engineering, Johns Hopkins University, 126 Latrobe Hall, 3400 North Charles Street, Baltimore, Maryland 21218 S. Bielamowicz Division of Otolaryngology, George Washington University, Washington, DC 20052 (Received 23 August 2010; revised 21 September 2011; accepted 23 April 2011) An immersed-boundary method based flow solver coupled with a finite-element solid dynamics solver is employed in order to conduct direct-numerical simulations of phonatory dynamics in a three-dimensional model of the human larynx. The computed features of the glottal flow including mean and peak flow rates, and the open and skewness quotients are found to be within the normal physiological range. The flow-induced vibration pattern shows the classical “convergent-divergent” glottal shape, and the vibration amplitude is also found to be typical for human phonation. The vocal fold motion is analyzed through the method of empirical eigenfunctions and this analysis indicates a 1:1 modal entrainment between the “adduction-abduction” mode and the “mucosal wave” mode. The glottal jet is found to exhibit noticeable cycle-to-cycle asymmetric deflections and the mechanism underlying this phenomenon is examined. C 2011 Acoustical Society of America. [DOI: 10.1121/1.3592216] V PACS number(s): 43.70.Bk, 43.70.Aj, 43.70.Gr [DAB] I. INTRODUCTION Phonation results from a complex nonlinear biomechanical coupling between glottal aerodynamics and vocal fold tissue. Computational modeling of phonation usually involves two key constituents, the model for the glottal aerodynamics and the model for vocal-fold (VF) deformation and vibration. In the past, various computational approaches have been proposed and employed in these two arenas. Models for the vocal fold dynamics range from the simplest “one-mass” (Flanagan and Landgraf, 1968), two-mass (Ishizaka and Flanagan, 1972) and multi-mass models (Story and Titze, 1995; Titze, 1973, 1974; Yang et al., 2010) to more accurate “continuum-based model (Rosa et al., 2003; Tao et al., 2007; Zheng et al., 2009). Similarly models for the glottal aerodynamics range from the simplest Bernoulli equation based model (Berg et al., 1957; Flanagan and Landgraf, 1968; Ishizaka and Flanagan, 1972; Story and Titze, 1995), to modified inviscid Euler (Lamar et al., 2003) and steady Navier-Stokes equation (Rosa et al., 2003). Unsteady Navier-Stokes equation based models which offer the highest level of fidelity (Tao et al., 2007; Duncan et al., 2006; Zheng et al., 2009) are also being used routinely in this arena. Increasing modeling fidelity enables us to capture additional features of phonation and provides better insights into the underlying biophysics. For example, the initial “onemass model” coupled with the modified Bernoulli equation was able to produce sustained vibrations. Subsequent development of the “two-mass model” enabled the reproduction a) Author to whom correspondence should be addressed. Electronic mail: mittal@jhu.edu 404 J. Acoust. Soc. Am. 130 (1), July 2011 Pages: 404–415 of the mucosal wave and a reasonable glottal wave form. Similarly, while Bernoulli equation based flow models may capture some basic features of the glottal pressure variation, only Navier-Stokes based models generate features such as asymmetric glottal jet deflection (Tao et al., 2007, Zheng et al., 2009) With regards to Navier-Stokes flow models, it should be noted that the vast majority of the past studies have employed two-dimensional models (Zhao et al., 2002; Zhang et al., 2002; Vries et al., 2003; Tao et al., 2007; Duncan et al., 2006; Luo et al., 2008; Zheng et al., 2009; Link et al., 2009). While these two-dimensional models, which represent the flow and VF vibration across a nominal anteriorposterior glottal midline, have provided good insights into phonation, there are several limitations associated with them. First, high speed imaging studies have shown that the flow-induced vibration pattern is three dimensional with a  strong anterior-posterior variation (Svec and Schutte, 1996). Second, the resulting glottal flow also has a highly threedimensional structure. Recent experimental studies have infact shown that the glottal flow is highly three-dimensional even within an idealized rectangular duct model of the larynx (Triep and Brücker, 2010). The three-dimensional jet structure is a direct result of vortex stretching and tilting which are also responsible for the eventual transition to turbulence of the glottal jet Finally, the sub and supraglottal lumen in the human larynx has significant variations in the anterior-posterior (z) direction which likely produce additional three-dimensional effects in the glottal flow. Threedimensionality in the flow structure could modify pressure losses in the glottis as well as the voice generation and its spectrum. Finally, three-dimensionality in the vibratory pattern has implication for laryngeal pathologies such as vocalfold nodules and polyps (Gunter, 2004). 0001-4966/2011/130(1)/404/12/$30.00 C 2011 Acoustical Society of America V While the search for biophysical insights drive most of the development of computational models in biomechanics, rapid advancements in computational hardware and software capabilities raise the possibility of conducting surgical planning based on such models. This prospect is particularly promising for phonosurgery due to the fact that phonation is primarily governed by biomechanical phenomena, which are amenable to model based prediction. Medialization laryngoplasty (Bielamowicz and Berke, 1995) is one surgical procedure, that is, particularly well suited for model based surgical planning. The procedure is commonly used to restore voice in patients with unilateral vocal-fold paralysis/ paresis and involves the insertion of a synthetic implant into the larynx via a window cut into the thyroid cartilage. A properly placed implant medializes the involved vocal fold thereby leading to sustained vocal-fold vibrations and restoration of the voice. The surgical outcome is however highly dependent on the shape, size and location of the implant, and the surgeon has few tools at hand that can predict the effect of different implant configurations on a given patient. A computational biomechanics tool that provides such a capability would significantly advance the surgical management of laryngeal pathologies and potentially improve surgical outcome and reduce revision rate. The level of fidelity required in such a surgical tool necessitates the use of sophisticated modeling approaches including continuum models for vocal fold mechanics and three-dimensional Navier-Stokes based glottal flow modeling. Furthermore, simulations would have to be performed in patient-specific three-dimensional geometries. A necessary precursor to the development of such a tool is to demonstrate the ability to conduct three-dimensional simulations in simple laryngeal geometries and to generate flow and vibration characteristics that can be verified against established experimental results. Such three-dimensional (3D) simulations are currently extremely challenging due to the computational expense involved; an adequately resolved three-dimensional simulation of phonatory dynamics can be hundred-times more expensive than a corresponding two-dimensional simulation. This has severely limited the deployment of 3D laryngeal computational models. To the best of our knowledge, only one three-dimensional flow-structure interaction study which coupled a Navier-Stokes flow model with a continuum vocal fold model has been documented in phonation research to date. In this study of Rosa et al., (2003), a 3D continuum, transversally isotropic, viscoelastic vocal fold model was coupled with a 3D, steady Navier-Stokes equation based flow model, and the simulation was performed in a cylindrical pipe which included false vocal folds. The simulation did manage to capture self-sustained divergent and convergent vocal folds vibration. However no details regarding the flowfield were provided. This was possibly due to low grid resolution (2600 tetrahedrons for the airflow and 3000 tetrahedrons for the vocal folds tissue) employed in this study or due to the assumption of steady flow, which eliminated all intrinsic time-dependent features of the flow. The objective of the current study is to document the results of a three-dimensional, fully unsteady [as opposed J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 the quasi-steady approach such as that of Rosa et al., (2003)], biomechanical model of phonation. We employ a relatively simple geometrical model of the larynx, that is, similar to some past experimental studies (Erath and Plesniak, 2006a,b, 2010; Becker et al., 2009), and the focus of the current study is on a critical evaluation of the computed results and verification against available experimental and numerical data. The simulations also allow us to gain some new insights into glottal jet and vocal fold dynamics and these are also discussed here. II. METHOD The computational method for flow-structure interaction, that is, employed in this study is briefly described here. This includes an immersed boundary method (IBM) method for the glottal airflow modeling, a finite-element method for the vocal fold mechanics and an efficient coupling procedure for the two. This numerical algorithm has been discussed in greater detail in Zheng (2009). A. Glottal airflow modeling The equations governing the glottal flow are the unsteady incompressible Navier–Stokes equations @vi ¼ 0; @xi @vi @vi vj 1 @p @ 2 vi þ ; þt ¼ @t @xj @xj q @xi @xj (1) where vi are velocity components, p is pressure, and q and t are flow density and kinematic viscosity. The equations are discretized in space using a cell-centered, non-staggered arrangement of the primitive variables vi and p. A second-order fractional step method is used to integrate the equations in time. The advection term is linearized using a second-order Adams-Bashforth scheme. Also, an implicit Crank-Nicolson scheme is employed to discretize the diffusion term and to eliminate the viscous stability constraint. A line-successive-overrelaxation (SOR) scheme is used to solve the advection-diffusion equation and an alternating-direction, geometric multi-grid with a line-SOR smoother is used to solve the pressure Poisson equation. The boundary conditions are implemented through a sharp-interface immersed-boundary method (IBM), which is locally and globally second-order accurate in space, and well suited for the simulation of complex and moving boundary problems on a Cartesian grid. The details of this flow solver as well as validation and verification studies of the methodology can be found in Mittal et al., (2008). B. Vocal fold tissue modeling The human vocal folds primarily consist of the cover, ligament, and vocalis (Hirano, 1977). From a biomechanical view-point, the vocal folds can be assumed to be made of a viscoelastic material which is transversally isotropic (Alipour et al., 2000). Since during normal phonation, vocal folds vibration induces a deformation, that is, small (Alipour et al., 2000), the vocal folds may be considered to behave as Zheng et al.: Three-dimensional laryngeal model 405 FIG. 1. (Color online) Anatomical model reconstruction of vocal folds and false vocal folds and the three layer inner structure of vocal folds. (a) A coronal view of human larynx from a high resolution CT scan. (b) An idealized shape three-dimensional vocal folds and false vocal folds. (c) Computational domain employed in the current study. a linear viscoelastic material. The constitutive law for this type material is as follows (Belytschko et al., 2000) rij ¼ Cijkl ekl þ Aijkl e_kl ; (2) where r is stress tensor, e is train tensor, e_ is strain rate tensor, and C and A are fourth-order tensor material constants. The dynamics of continuum vocal folds are governed by the Navier equation: rij;j þ qfi ¼ q€ ui ; (3) where i and j range from 1 to 3, fi is the body force component in i direction, q is vocal fold tissue density, ui is displacement in i direction, and u€i is the acceleration component in i direction. This equation is solved by the finite-element method (FEM). Integration of the weak form of the above equation with subsequent application of the divergence theorem results in the following discrete equation: Mab U€b þ Cab U_ b þ Kab Ub ¼ Fa ; J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 C. Flow-Structure Interaction The large (1000) solid-fluid density ratio in phonation allows effective use of explicit coupling (Zheng, 2009) which is adopted in the current model. The flow and solid solvers communicate through the Lagrangian interface which is represented by the triangular surface mesh which is shown in Fig. 1(b). In the explicit-coupling method, the flow is marched by one step with the existing deformed shape and velocities of the solid as the boundary conditions. The aerodynamic forces on the VF are then computed at the vocal fold surface through an interpolation scheme on the flow grid. Finally, the equation for the solid is marched by one time step with the updated surface traction, and the deformation and velocities on the solid grid are transferred to the vocal fold surface, so that the fluid/solid interface can be updated. (4) where M represents the mass matrix, C represents the damping matrix, K represents the stiffness matrix, F represents finite element nodal force, and U represents the nodal displacement. Equation (4) is a second-order ordinary differential equation in time and is discretized using the second-order Newmark scheme (Belytschko et al., 2000) wherein a banded LU decomposition is used to solve the system of algebraic equations. The Cuthill-McKee (Cuthill and Mckee, 1969) and Gibbs-Poole-Stockmeyer methods (Gibbs et al., 406 1976) are used to re-index the nodes in order to create banded matrixes. The details of this solid solver can be found in Zheng (2009). III. SIMULATION SETUP In the current study, a three dimensional idealized shaped vocal fold model is employed which is 1.5 cm long in the anterior-posterior (z-) direction and is invariant in shape in the longitudinal direction [shown in Fig. 1(b)]. The vocal fold profile on the cross section is nominally based on a coronal view [shown in Fig. 1(a)] of a high resolution laryngeal CT scan of a 30 year old normal male subject. Due to the slight asymmetry of the subject’s larynx, only the right part is used and the left part is created by mirroring the right Zheng et al.: Three-dimensional laryngeal model TABLE I. Material properties of the three-layers of the vocal folds. Property Layer Cover Ligament Body qðg=cm3 Þ Ep ðkPaÞ vp Epz ðkPaÞ vpz Gpz ðkPaÞ gðPoiseÞ 1.043 1.043 1.043 2.041 3.306 3.990 0.9 0.9 0.9 20 000 33 000 40 000 0.0 0.0 0.0 10 40 20 1 1 1 part about the centerline. Both the true and false vocal folds are represented and the resulting dimensions of the larynx are as follows: The true vocal folds (TVFs) are 1 cm long and extend 0.99 cm toward the glottal midline. The false vocal folds (FVFs) are 2.3 cm long and extend 0.67 cm toward the supraglottal space. The ventricles are about 0.56 cm wide at their widest location and the false vocal fold gap is 0.667cm. These dimensions and shape are the same as those used in our previous two-dimensional study (Zheng et al., 2009) and shown to be in the range typical for human adults. The internal anatomical structure of the vocal folds is made up of three layers (Hirano, 1977): cover, ligament and vocalis. The precise configuration of the layers is adopted from Luo et al., (2008) and each layer is assumed to be invariant in shape in anterior-posterior direction [as shown in Fig. 1(b)]. The material properties for each layer are taken from Alipour et al., (2000) and are shown in the Table I where q is the tissue density, Ep is the transversal Young’s Modulus, cp is the in-plane transversal Poisson’s ratio, Epz is the longitudinal Young’s Modulus, cpz is the longitudinal Poisson’s ratio, Gpz is the longitudinal shear modulus. It should be noted that due to the in-plane motion assumption, the longitudinal Young’s modulus Epz is not given in Alipour et al., (2000). In the current study, the longitudinal Young’s Modulus Epz is set to equal to 104 times the transversal Young’s modulus Ep and this automatically enforces the satisfaction of the in-plane motion constraint (Cook et al., 2008). Since the false vocals hardly move during normal phonation, they are treated as rigid bodies. A 12 cm  2.0 cm  1.5 cm straight rectangular duct is used to represent the human larynx lumen [shown in Fig. 1(c)] and this is similar to many past laboratory experiments (Triep et al., 2005; Erath and Plesniak, 2006a,b; Becker et al., 2009). The true vocal folds are located from x ¼ 2.16 cm to 3.16cm in the anterior-posterior (z-) direction, which results in a location of the glottal exit exactly at x ¼ 3 cm. A zero gauge pressure is prescribed at the outlet and the gauge pressure at the inlet is set to a typical value of 1 kPa (Titze, 1994). No-slip, no-penetration flow velocity conditions are applied on all the walls. For the solid vocal folds, zero displacement boundary conditions are prescribed on the anterior (z ¼ 0 cm) and posterior (z ¼ 1.5 cm) ends as well as the lateral walls (y ¼ 0 cm for the lower vocal fold and y ¼ 2 cm for the upper vocal fold). The traction boundary conditions (computed aerodynamic loads) are applied on the vocal fold surfaces. A high-resolution non-uniform 256  128  64 Cartesian grid is employed for the flow solver and a 58 427 tetrahedral element grid is employed for the solid solver. The grid is based on our experience with 2D simulations of similar configurations (Zheng et al., 2009) and provides the highest resolution in the intraglottal and near supraglottal region where high flow gradients are expected. A small time step of 1:75  103 ms is employed to provide good temporal resolution as well as to satisfy the CFL stability constraint. This simulation is carried out on an IBM iDataplex parallel computer using 128 Intel 3.0 GHz Quad-Core processors. The wall-clock time for the simulation is about 720 h and the overall computational expense of one simulation is over 90 000 CPU hours. IV. RESULT AND DISCUSSION A. Glottal volume velocity waveform Figure 2(a) shows the time-history of the glottal flow rate which indicates a short transient period that quickly FIG. 2. (a) Computed time-history of the glottal flow rate. (b) Schematic of the typical glottal airflow waveform (adopted from Titze, 1994). (c) Glottal airflow wave form computed in current study for four different cycles. J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 Zheng et al.: Three-dimensional laryngeal model 407 develops into a fairly stationary limit-cycle vibration. In the transient stage, a strong sweeping vocal fold motion is observed but is found to damp out rapidly. The vocal fold vibration pattern in the transient stage as well as in the stationary state will be discussed in details in next section. It should be also noted that during the stationary state limitcycle vibration, the peak glottal flow rate shows some cycleto-cycle variability. This irregular variation has been also reported in the stiffer vocal fold cases of Rosa et al., (2003) and could be due to the nonlinear effects resulting from the complex flow-structure interaction. The typical mean glottal flow rate measured in the experiments ranges from 110 ml/s to 220 ml/s (Zemlin, 1998) and the peak glottal flow rate ranges 200 ml/s to 580 ml/s (Hertegard et al., 1992). In the current study, the mean glottal flow rate and phase averaged peak glottal flow rate over the 15 stationary vibration cycles are 121 ml/s and 298 ml/s, respectively, which are within these typical ranges. The fundamental phonation frequency computed using the glottal flow rate is 243 Hz, which leads to a 4.2 ms vibration period. It should be noted that this frequency is relatively high, but still within the human phonation frequency range from 60 Hz to 260 Hz (Zemlin, 1998). The relatively high phonation frequency here is a consequence of the large longitudinal Young’s modulus and is inline with Rosa et al., (2003) who noted an increase in phonation frequency from 160 Hz to 190 Hz when the longitudinal Young’s modulus was increased from 20 kPa to 200 kPa. The longitudinal stiffness incorporates the effect of the muscle activation and modeling this effect remains a challenge for these models. According to the source-filter theory (Fant, 1960), the sound source can be represented by the time-varying glottal airflow waveform, that is, filtered by the vocal tract. A typical glottal airflow waveform is shown in Fig. 2(b), which can be characterized by several parameters that are important determinants of sound quality (Titze, 1994). As shown in Fig. 2(b), the vibration period is defined as T, Tp is time of increasing flow, Tn is time of decreasing flow, and T0 is the duration of flow. Two dimensionless parameters s0 and ss are defined from these quantities, where s0 is the open quotient defined as T0 =T and ss is skewing quotient defined as Tp =Tn . Physiologically, the open quotient Q0 can be adjusted by the glottal width and typical values of Q0 range from 0.4 to 0.7 (Titze, 1994). A value of Q0 which is lower than 0.4 indicates a “pressed” sound and a value higher than 0.7 indicates a “breathy” sound (Titze, 1994). For the current study, the average Q0 is 0.62, which is within the normal range. It should be noted that a minimum glottal gap of two gridpoints, which corresponds physically to a gap of 0.16 mm (or about 0.8% of the channel width) is maintained throughout the simulation due to the requirements of the flow solver (Mittal et al., 2008). This gap produces a small amount of flow leakage equivalent to about 4% of the peak flux, even during full closure and likely increases the open quotient slightly. It is interesting that incomplete closure can oftentimes be observed during normal phonation in healthy larynges (Titze, 1994) although leakage flux rates have not been measured. The small amount of leakage is not expected to have any significant effect on the glottal jet dynamics. The typical glottal flow wave form exhibits a slow rise followed by a rapid fall. This skewness results from the delay of the buildup of the flow with respect to the movement of the vocal folds due to the inertia of the air column and the sudden collapse of vocal folds (Titze, 1994). The skewness quotient shows a wide range from 1.1 to 3.4 (Ishizaka and Flanagan, 1972; Lamar et al., 2003; Duncan et al., 2006; Luo et al., 2008) and the current value of 1.32 is well within this observed range. B. Vocal fold vibration The vocal fold vibration results from the nonlinear interaction of glottal flow and vocal fold tissue. In past studies, the vocal fold vibration was usually captured from a superior aspect using in vivo endoscopic recordings (Baer, 1975; Krausert et al., 2010). Due to the optical inaccessibility of the inferior part of the vocal folds, the three-dimensional details of the vocal motion fold motion/deformation remain unclear. During phonation, the vocal folds are believed to exhibit a characteristic “converging-diverging” vibratory pattern. This type of motion has also been shown to result from “mode entrainment” (Zemlin, 1998) of two distinct vibratory modes; the so called “adduction and abduction” mode and the “mucosal wave” mode. Figure 3 shows the vocal-fold vibration pattern captured in this study, which includes the superior views as well as the center-plane configurations of the vocal folds at six FIG. 3. (Color online) Vocal fold shape at six different time-instants within a vibration cycle. For each time instance, the figure on the left and right show the superior view and the medial vocal fold profiles, respectively. (a) 0.06562s (12%T), (b) 0.006685s (42%T), (c) 0.06755s (58%T), (d) 0.06808s (70%T), (e) 0.06860s (83%T), (f) 0.06913s (95%T), Note that T ¼ 0.0042s is the vibration period and this particular cycle starts at 0.00651s. The corresponding glottal flux rates at these time instants can be found in Fig. 2(c). 408 J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 Zheng et al.: Three-dimensional laryngeal model FIG. 4. (Color online) The most energetic three empirical eigen functions of vocal fold motions. The first two rows are perspective views of the empirical modes at two extreme phases and the third row is the corresponding medial profiles at the center-plane (z ¼ 0.75) (a) mode-1, (b) mode-2: and (c), mode-3. different time instants within one vibration cycle. The “converging-diverging” vibration pattern can be clearly observed from the series of figures. As shown in Fig. 3(b) and 3(c), the vocal folds start to open from beneath at the beginning of the cycle and the glottis presents a convergent shape. In Fig. 3(d), the glottis reaches its maximum opening and starts to close first from the inferior portion followed by the superior [shown in Fig. 3(d) and 3(f)] portion. The propagation of the mucosal wave can be clearly observed and the largest vibration is found to occur at the superior tip of the vocal folds. It is interesting to note that the same patterns were reported by Berry et al., (2001) using high-speed digital imaging in an excised canine hemilarynx. Thus a physically viable and qualitatively verifiable vibration pattern is successfully captured in the current study. The computed vocal fold deformation allows for a detailed analysis of the vibratory pattern. Here we employ proper orthogonal decomposition (POD) analysis to extract the most energetic empirical eigen-modes. First, the time-mean vocal fold profiles and displacement time variations are computed over 15 vibration cycles from t ¼ 0.02 s to 0.0814 s, which contains 351 equally spaced snapshots of vocal fold profiles. Let uk ðxÞ represents the kth snapshot of the displacement field where 1  k  N and N ¼ 351 is the total number of snapshots in the current study. The correlation matrix S associated with the ith and jth snapshots can be computed as follows (Lewin and Haj-Hariri, 2005) where U is the matrix of snapshots fu1 ; u2 …uN g, and km and qm are the mth eigenvalue and eigenvector, respectively, of the correlation matrix S. The modal coefficient can subsequently be obtained by ym ðkÞ ¼< uk ; um > : (7) The eigenvalue represents the amount of energy associated with the empirical eigenmodes and Fig. 4 shows the three most energetic modes, which together, capture 93% of the total energy. The first empirical mode corresponds to an adduction-abduction motion (medial-lateral vibration) and captures 50% of the total energy. The second mode is a convergent-divergent motion and contains 25% of total energy. The third mode is a “sweeping motion” (inferior-superior vibration) and captures 18% of total energy. The residual 7% energy is distributed into the rest of the modes. The modal coefficient time history is shown in Fig. 5. Through power spectrum analysis of these modal coefficients, we find that mode-1 and mode-2 vibrate at the phonation frequency of 243 Hz. This is a clear manifestation of 1:1 mode entrainment which has been found in past studies (Berry et al., 1994; Berry et al., 2001; Zhang et al., 2007). Berry et al., (2001) have also extracted the two-dimensional N Si;j ¼< ui ; uj > pt 1 X ui ðxl Þ  uj ðxl Þ; Npt l¼1 (5) where Npt are the total number of surface points. The empirical eigen-modes um can then be computed by pffiffiffiffiffiffi um ¼ 1= km Uqm ; J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 (6) FIG. 5. (Color online) Modal coefficient time history of the three most energetic empirical eigen-modes: Mode 1: solid line, Mode 2: dashed line, Mode 3: dashed-dotted. Zheng et al.: Three-dimensional laryngeal model 409 FIG. 6. (Color online) Glottal gap-width time history at five different longitudinal planes: Z ¼ 0.75: solid line, Z ¼ 0.6: dashed line, Z ¼ 0.9: dasheddotted, Z ¼ 0.375: dashed line, Z ¼ 1.125: dashed-dotted. empirical modes for a hemilarnxy. In their study, 1:1 modal entrainment was captured and a 3:1 modal energy ratio shown between the “adduct-abduct” mode and the “mucosal wave” mode. It should be noted that the modal energy ratio in the current study is 2:1 which is different from Berry’s study. However, Zhang et al. (2007) have shown that this modal energy ratio also depends on several parameters, such as structural damping and jet velocity. Thus the difference in energy partitioning should be expected from the studies with different parameters. Mode-3 in the current study is found to have two distinct frequencies of 147 Hz and 243 Hz although the 147 Hz frequency has a much larger energy density in the spectrum. The current three-dimensional study also allows us to access the anterior-posterior vocal fold motion. As shown by the superior views of vocal folds in Fig. 3, the maximum vibration amplitude always occurs at the center plane (z ¼ 0.75 cm) and decreases monotonically toward the two ends. Figure 6 shows the glottal gap width time history measured at five different anterior-posterior planes. It is found that there is no phase difference along the longitudinal direction and furthermore, the vibration is always symmetric about the center plane. A similar vibration pattern has been captured in the recent experimental study of Becker et al., (2009), which also used an idealized vocal fold model. This longitudinal symmetry is likely a result of the glottal shape invariance along the longitudinal direction (Berry et al., 1994). We also note that the maximum glottal gap width is about 0.1 cm, which is inline with previous two-dimensional numerical studies (Tao et al., 2007; Luo et al., 2008; Zheng et al., 2009). C. Glottal jet dynamics The glottal flow exhibits a variety of interesting phenomena such as jet instability, flow transition, turbulence and asymmetric jet deflection etc. all of which have implication for phonation and the quality of voice. Direct numerical simulation of the glottal flow provides us with the opportunity to perform a detailed analysis of the spatial and temporal features of the flow. Figure 7 shows the iso-surface of FIG. 7. (Color online) The iso-surface of enstrophy at four different time-instants within four different vibration cycles. (a) the 11th cycle with a downward jet deflection, (b) the 12th cycle with a upward jet deflection, (c) the 17th cycle with a centered jet, (d) the 21st cycle with split jet deflection. 410 J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 Zheng et al.: Three-dimensional laryngeal model enstrophy at four time-instants during four different vibration cycles. The jet is nearly straight when the flow is initially pushed out of the glottis (the first row of Fig. 7). As the glottis continues to open, the length of the jet column increases and during this stage, the jet is found to deflect from the centerline. The jet deflection is observed to vary significantly from cycle to cycle; For example, the jet exhibits a downward deflection in the 11th cycle [shown in Fig. 7(a)], an upward deflection in 12th cycle [shown Fig. 7(b)], and is nearly undeflected in the 17th cycle [shown in Fig. 7(c)]. The jet deflection also sometimes exhibits an anteriorposterior variation; for example, in the 21st cycle [shown in Fig. 7(d)], the jet is deflected downward on one side and deflected upward on the other side. Glottal jet deflection has been described in several recent experimental studies (Erath and Plesniak, 2006a,b, 2010; Triep et al., 2005). The underlying physical mechanism is usually associated with the so called “Coanda effect” (Erath and Plesniak, 2006; Tao et al., 2007; Becker et al., 2009) which is intimately linked to the divergent shape of the glottis. Recently, we have conducted an extensive study of asymmetric glottal jet deflection using a two-dimensional laryngeal model (Zheng et al., 2011). The key conclusion of the study was that intrinsic instabilities of the supraglottal flow produce asymmetric flow patterns that propagate upstream to the glottis and eventually produce large-scale glottal jet deflections. However, a limitation of this study was the use of a two-dimensional model that eliminated all three-dimensional effects. Thus the current three-dimensional model allows us to evaluate the validity of the previous two-dimensional results In order to quantify and further analyze the jet deflection, we focus on the jet immediately downstream of the glottis (x ¼ 3.1cm). Assuming that the center of jet coincides with the location of the maximum velocity magnitude, ~ max , the jet deflection angle (denoted by h) at any longitujVj dinal plane can be represented by the tangent to the streamline at this location. Thus,   V  (8) max; h ¼ tan1 U jVj ~ where U and V are the streamwise (x-direction) and mediallateral (y-direction) velocity components, respectively. We also compute this angle only during the phase of the glottal cycle when the glottis is open and there is measurable flow. For the rest of the cycle, there is almost no flow and the deflection angle is set to be zero. Figure 8(a) shows the time history of the jet deflection computed at the jet centerplane (z ¼ 0.75 cm) for all of the computed cycles. The glottal gap width is superimposed in the figure for ease of discussion. A number of insightful observations can be made from this plot. First, it is noted that the jet deflection during the first four cycles is very small but starting from the fifth cycle, a clear pattern of the significant jet deflection is established. This behavior is inline with our observation for a 2D laryngeal model (Zheng et al., 2011) and supports our assertion that the glottal jet deflection is driven by remnant J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 FIG. 8. (a) Time history of the measured jet deflection angle at center plane and downstream of the glottis solid line: jet deflection angle; dashed line: glottal gap. (b) Comparison of jet deflection at three anterior-posterior (Z direction) locations. asymmetric vortex structures which develop downstream in previous cycles and take a few cycles to propagate back to the glottis. Beyond the fifth cycle, we note that during the glottal opening (flow acceleration) phase, the deflection angle increases and quickly reaches a maximum which varies from 10 to 25 . The deflection then starts to decrease and reduces almost back to 0 at the maximum glottal opening phase. Subsequently, during the closing phase (flow deceleration), the jet deflection starts to oscillate randomly about 0 with a small amplitude. The initial large deflection followed by a straightening of the jet at maximum opening is also similar to that observed with the 2D model. This further supports the notion that the glottal jet deflection is driven by downstream flow structures (Zheng et al., 2011) since deflections are expected to be smallest when the jet velocity is large enough to overwhelm the disturbance due to the remnant vortices. Interestingly, the relatively small jet deflection during the flow deceleration phase observed here is in disagreement with the 2D simulations (Zheng et al., 2011) where the largest deflections are usually observed during this phase. In the two-dimensional flow model, the jet is most unstable during the deceleration stage (Zheng et al., 2011) and this coupled with the relatively low jet velocity makes the jet most susceptible to large-scale deflections at this stage. In contrast, in the three-dimensional flow model, as the glottis opens and the jet emanates from the glottis, three-dimensional instabilities are initiated and jet immediately downstream of the glottis becomes increasingly threedimensional. This three-dimensionality results in reduced anterior-posterior (z-direction) coherence in the jet and makes it more resistant to large-scale jet deflection. Note Zheng et al.: Three-dimensional laryngeal model 411 FIG. 9. (Color online) Iso-surface of magnitude of span wise vorticity at four different time-instants within one cycle. that although this behavior is contrary to what is seen in the 2D simulations, it nevertheless further contradicts the notion that jet deflection is driven by the Coanda effect. This is because the largest deflections occur during the glottal opening phase when the glottis has a convergent shape, and almost no deflections occur when the glottis has a divergent shape. It should be noted that the importance of downstream flow structures in driving jet deflection has also been suggested in the recent studies of Neubauer et al., (2007) and Triep and Brücker (2010). Figure 8(b) shows the jet deflection at three different anterior-posterior (z-) locations for four consecutive cycles. This plot further demonstrates the effect of three-dimensionality on jet deflection. Of the four cycles shown here, we note that the jet deflects in the same direction at all three points in the second and third cycles. However for the first and fourth cycles, the jet does not deflect in the same direction at all the locations, indicating significant three-dimensionality in the deflection of the jet. Figure 9 shows the top view of anterior-posterior vorticity (xz ) at four different time-instants within one cycle, and the top vocal and false folds are removed to more clearly show the flow features. The jet flow in the longitudinal direction exhibits a time-varying structure which is strongly three dimensional. When the jet is initially just pushed out of the glottis, the flow is almost uniform along the longitudinal direction [Fig. 9(a)]. However, as more flow is pushed out of glottis and the jet column increases in length, the flow starts to become three dimensional. Furthermore, the anterior-posterior (z-direction) jet-width decreases along the streamwise (x-) direction and forms a “V”-shaped pattern [Fig. 9(b) and 9(c)]. Triep and Brücker (2010) have documented the same type of flow pattern using PIV and have concluded that this vena contracta leads to axis switching of the jet This is 412 J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 FIG. 10. (a) Five locations along the center line in x-y plane chosen to perform the flow spectrum analysis. (b) Span-averaged stream-wise velocity spectrum at the five different locations: (1) x ¼ 3.025 cm, (2) x ¼ 3.175 cm, (3) x ¼ 3.325 cm, (4) x ¼ 3.575 cm, and (5) x ¼ 3.975 cm, respectively. Dashed line corresponds to k5=3 . therefore another computed flow feature that can be verified against experiments. The sound production mechanism in phonation is directly related to the features of the postglottal flow including the transition to turbulence. In order to further understand the transition to turbulence in the glottal jet, we examine the frequency spectra of the velocity in the jet core. We have computed the frequency spectra at the five locations along the center line [shown in Fig. 10(a)]. The spectra are also averaged along the anterior-posterior direction (z-) in the core of the glottal jet so plotted in the spectra is a line with a slope of 5/3 which represents the inertial subrange (Tennekes and Lumley, 1972). The plots in Fig. 10(b) show the clearest presence of an inertial subrange at location 4 indicating that the transition to turbulence occurs between locations 3 and 4. It should be noted that the peak jet Reynolds number [defined following Triep and Brücker (2010) as Re ¼ UWmax =, where U the peak flux averaged velocity, Wmax the maximum glottal gap width and  the kinematic viscosity of air] is about 210, which is one-fifth of the typical Reynolds number for an adult human. The lower Reynolds number is necessitated by the resolution requirements and associated computational expense of these direct-numerical simulations. However, despite this lower Reynolds number, the current simulations allow us to capture the phenomenon of transition to turbulence in the glottal jet. It is likely that at the higher physiological Reynolds numbers, the transition will occur closer to the glottis and lead to a higher turbulence intensity in the glottal jet. Zheng et al.: Three-dimensional laryngeal model FIG. 11. (Color online) Phase-averaged pressure distribution measured at two different directions (a) x is varied from 0 to12 cm for y ¼ 1.0 cm and z ¼ 0.75cm, (b) z is varied from 0 to 1.5 cm for x ¼ 2.975 cm and y ¼ 1.0 cm, (c) Isosurface of pressure (p ¼ -0.1248 kPa) at t ¼ 0.0441 s. (d) Four locations along the center line in the center plane (z ¼ 0.75 cm) where pressure time-history is extracted, (e) Pressure time history at four different locations and corresponding glottal flow rate. D. Glottal pressure distribution The variation of intraglottal and supraglottal pressure is important from the standpoint of understanding self-sustained vocal fold vibrations as well for developing/parameterizing lower fidelity models, and the current simulations allow us to examine the effect of three-dimensionality on the pressure field. Figure 11(a) shows the phase averaged pressure distribution along the streamwise (x-) direction at the center of the channel (y ¼ 1.0 cm; z ¼ 0.75 cm). At the flow acceleration phase (0–0.4T), there is a favorable pressure gradient in the supraglottal region which helps convect the flow away from the glottis. In contrast, during the deceleration phase (0.4T–0.55T), a strong adverse pressure gradient develops in the supraglottal region which enhances local flow instabilities. These distributions have some similarity to past studies done with static vocal fold models (Scherer et al., 2001, 2002) which would indicate that dynamic effects due to VF vibration might have a relatively minor effect on the transglottal pressure distribution. Figure 11(b) shows the phase averaged pressure distribution along the anterior-posterior (z-) direction at a location just inside the superior glottal opening at x ¼ 2.975 cm and y ¼ 1.0 cm. It can be observed that the pressure varies significantly along the anterior-posterior direction which is a clear manifestation of 3D effects on the pressure. Further evidence of the effect of three-dimensionality on the pressure can be seen in Fig. 11(c) which shows an isosurface of pressure, corresponding to p ¼ 0.01248 kPa at t ¼ 0.0441s when the J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 two vocal folds are fully open. The pressure distribution at this time-instant is highly three-dimensional and the topology corresponds well with the corresponding vortex topology at this instant show in Fig. 7. The pressure wave forms (time histories) are also shown at four different locations along the center line in the center plane (z ¼ 0.75cm), which are located at the entrances and exits of glottis and ventricles, respectively [shown in Fig. 11(d)]. As shown in Fig. 11(e), at the beginning of flow acceleration phase, the pressure decreases first and reaches a minimum value. It subsequently starts to increase during the flow deceleration phase, reaching a peak value before decreasing again. A slight phase difference is observed between the upstream location (x ¼ 2.57 cm) and the downstream locations (x ¼ 3.08, 3.48, 4.23 cm) which is expected due to the acceleration of the flow through the glottis. E. Conclusions A direct numerical simulation of flow in a 3D laryngeal model that includes a fully coupled, continuum based structural vocal fold model, has been carried out to study the glottal flow and vocal fold vibration during phonation. The resulting glottal flow rate and wave form are found to be within the normal physiological range. The simulation captures the convergent-divergent vibratory pattern of the vocal folds, that is, characteristic of normal phonation (Titze, 1994). POD analysis shows that the abduction-adduction mode is the most energetic mode. This mode vibrates at the Zheng et al.: Three-dimensional laryngeal model 413 same frequency as the mucosal-wave type mode which indicates that 1:1 modal entrainment (Berry et al., 1994; Berry et al., 2001) has been successfully captured. A detailed analysis of the glottal jet supports the conclusion of our previous two-dimensional study (Zheng et al., 2011) in that the stochastic cycle-to-cycle jet deflections are a result of downstream flow structures and not some “Coanda effect” associated with the intra-glottal flow. The inclusion of threedimensionality, however, changes some aspects of the jet deflection including a significant reduction in jet deflection during the glottal closure phase. The simulated flow also exhibits a longitudinal vena contracta, which matches well with a recent experimental study (Triep and Brücker, 2010), and the simulation also captures the transition to turbulence in the glottal jet. The three-dimensional model developed here therefore produces many of the characteristic features of the glottal flow and vocal fold dynamics that have been documented in past studies. The simulations also provide new insights regarding jet deflection and transition. It should be noted that while these simulations are computationally expensive at the present time, trends in high-performance computing (Dongarra, 2004) suggest that such simulations will run hundreds (if not thousands) of times faster by the end of the decade. Ongoing extension of the current work is focused on adding more geometrical details to the laryngeal model with the ultimate objective of simulating phonation in a CT-based, patient-specific laryngeal model. ACKNOWLEDGMENTS The project described was supported by Grant Number ROIDC007125 from the National Institute on Deafness and Other Communication Disorders (NIDCD). The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIDCD or the NIH. This research was also supported in part by the National Science Foundation through TeraGrid resources provided by NICS under grant number TG-CTS100002. Alipour, F., Berry, D. A., and Titze, I. R. (2000). “A finite-element model of vocal-fold vibration,” J. Acoust. Soc. Am., 108(6), 3003–3012. Baer, T. (1975). “Investigation of phonation using excised larynxes,” Ph.D. dissertation, MA Institutes of Technology, MA. Becker, S., Kniesburges, S., Muller, S., Delgado, A., Link, G., Kalternbacher, M., and Dollinger, M. (2009). “Flow-structure-acoustic interaction in a human voice model,” J. Acoust. Soc. Am. 125(3), 1351–1361. Belytschko, T., Liu, W., and Moran, B. (2000). Nolinear Finite Elements for Continua and Structures (John Wiley, New York), p. 650. Berg, J. V. D., Zantema, J. T., and Doornenbal, P. J. (1957). “On the air resistance and the Bernoulli effect of the human larynx,” J. Acoust. Soc. Am. 29(5), 626–631. Berry, D. A., Montequin, D. W., and Tayama, N. (2001). “High-speed digital imaging of the medial surface of the vocal folds,” J. Acoust. Soc. Am. 110(5), 2539–2547. Berry, D. A., Herzel, H., Titze, I. G., and Krischer, K. (1994). “Interpretation of biomechanical simulations of normal and chaotic vocal fold oscillations with empirical eigenfunctions,” J. Acoust. Soc. Am. 95(6), 3595–3604. Bielamowicz, S., and Berke, G. B. (1995). “An improved method of medialization laryngoplasty using a three-sided thyroplasty window,” Laryngoscope 105(5), 537–539. Cook, D. D., Nauman, E., and Mongeau, L. (2008). “Reducing the number of vocal fold mechanical tissue properties: Evaluation of the incompressi414 J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 bility and planar displacement assumptions,” J. Acoust. Soc. Am. 124(6), 3888–3896. Cuthill, E., and McKee. J. (1969). “Reducing the bandwidth of sparse symmetric matrices,” in Proceedings of the 24th National Conference of the ACM (Association for Computing Machinery, New York), pp. 157–172. Dongarra, J. (2004). “Trends in high performance computing,” Comput. J. 47(4), 399–403. Duncan, C., Zhai, G., and Scherer, R. (2006). “Modeling coupled aerodynamics and vocal fold dynamics using immersed boundary methods,” J. Acoust. Soc. Am. 120(5), 2859–2871. Erath, B. D., and Plesniak, M. W. (2006a). “An investigation of bimodal jet trajectory in flow through scaled models of the human vocal tract,” Exp. Fluids. 40, 683–696. Erath, B. D., and Plesniak, M. W. (2006b). “The occurrence of the coanda effect in pulsatile flow through static models of the human vocal folds,” J. Acoust. Soc. Am. 120(2), 1000–1011. Erath, B. D., and Plesniak, M. W. (2010). “An investigation of asymmetric flow features in a scaled-up driven model of the human vocal folds,” Exp. Fluids 49(1), 131–146. Fant, G. (1960). Theory of Speech Production, (Mouton, The Hague), p. 328. Flanagan, J. L., and Landgraf, L. L. (1968). “Self-oscillating source for vocal-tract synthesizers,” IEEE Trans. Audio Electroacoust. 16(1), 57–64. Gibbs, N. E., Poole, W. G., and Stockmeyer, P. K. (1976). “A comparison of several bandwidth and profile reduction algorithms,” ACM Trans. Math. Software, 2(4), 322–330. Gunter, H. E. (2004). “Modeling mechanical stress as a factor in the etiology of benigh vocal fold lesions,” J. Biomech. 37(7), 1119–1124. Hertegard, S., Gauffin, J., and Karlsson, I. (1992). “Physiological correlates of the inversrsed filtered flow waveform” J. Voice. 6(3), 224–234. Hirano, M. (1977). “Structure and vibratory behavior of the vocal folds,” in Dynamic Aspect of Speech Production (University of Tokyo Press, Tokyo, Japan), pp. 13–27. Ishizaka, K., and Flanagan, J. L. (1972). “Synthesis of voiced sounds from a two-mass model of the vocal cords,” Bell Syst. Tech. J. 51(6), 1233– 1268. Krausert, C. R., Olszewski, A. E., Taylor, L. N., McMurray, J. S., Dailey, S. H., and Jiang, J. J., (2010). “Mucosal wave measurement and visualization techniques,” J. Voice, doi:10.1016/j.jvoice.2010.02.001. Lamar, M. D. Qi, Y., and Xin. J. (2003). “Modeling vocal fold motion with a hydrodynamic semicontinuum model,” J. Acoust. Soc. Am. 114(1), 455– 464. Lewin, G. C., and Haj-Hariri, H. (2005). “Reduced-order modeling of a heaving airfoil,” AIAA J. 43(2), 270–283. Link, G., Kaltenbacher, M., Breuer, M., and Döllinger, M. (2009). “A 2D finite-element scheme for fluid-solid-acoustic interactions and its application to human phonation,” Comput. Methods Appl. Mech. Eng. 198(41– 44), 3321–3334. Luo, H., Mittal, R., Bielamowize, S., Walsh, R., and Hahn, J. (2008). “An immersed-boundary method for flow-structure interaction in biological systems with applications to phonation,” J. Comput. Phys. 227, 9303– 9332. Mittal, R., Dong, H., Bozkuittas, M., Najiar, F. M., Vargas, A., and Loebbecke, A. (2008). “A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries,” J. Comput. Phys. 227(10), 4825–4852. Neubauer, J., Zhang, C., Miraghaie, R., and David, A. B. (2007). “Coherent structure of the near field flow in a self-oscillating physical model of the vocal folds,” J. Acoust. Soc. Am. 121(2), 1102–1118. Rosa, M. D. O., Pereira, J. C., Grellet, M., and Alwan, A. (2003). “A contribution to simulating a three-dimensional larynx model using the finite element method,” J. Acoust. Soc. Am. 114(5), 2893–2905. Scherer, R. C., Shinwari, D., Witt, K. J., Zhang, C., Kucinschi, R., and Afjeh, A. A. (2001). “Intraglottal pressure profiles for a aymmetric and oblique glottis with a divergence angle of 10 degrees,” J. Acoust. Soc. Am., 109(4), 1616–1630. Scherer, R. C., Shinwari, D., Witt, J. J., Zhang, C., Kucinschi, R., and Afieh, A.A. (2002). “Intraglottal pressure distributions for a symmetric and oblique glottis with a uniform duct (L),” J. Acoust. Soc. Am. 112(4), 1253–1256. Story, B. H., and Titze, I. R. (1995). “Voice simulation with a body-cover model of the vocal folds,” J. Acoust. Soc. Am. 97(2), 1249–1260.  Svec, J. G., and Schutte, H. K. (1996). “Videokymography: High-speed line scanning of vocal fold vibration,” J. Voice 10(2), 201–205. Zheng et al.: Three-dimensional laryngeal model Tao, C., Zhang, Y., Hottinger, D. G., and Jiang, J. J. (2007). “Asymmetric airflow and vibration induced by the Coanda effect in a symmetric model of the vocal fold,” J. Acoust. Soc. Am. 112(4), 2270–2278. Tennekes, H., and Lumley, J. L. (1972). A First Course in Turbulence (MIT Press, Cambridge), 266–276. Titze, I. R. (1973). “The human vocal cords: A mathematical model: Part I,” Phonetica, 28, 129–170. Titze, I. R. (1974). “The human vocal cords: A mathematical model: Part II,” Phonetica, 29, 1–21. Titze, I. R. (1994). Principles of Voice Production (Prentice-Hall, Englewood Cliffs), pp. 354. Triep, M., Brücker, C., and Schröder, W. (2005). “High-speed PIV measurements of the flow downstream of a dynamic mechanical model of the human vocal folds,” Exp. Fluids 39, 232–245. Triep, M., and Brücker, C. (2010). “Three dimensional nature of the glottal jet,” Acoust. Soc. Am. 127(3), 1537–1547. Vries, M. P., Hamburg, M. C., Schutte, H. K., and Verkerke, G. J. (2003). “Numerical simulations of self-sustained oscillation of a voice-producing element based on Navier-Stokes equations and the finite element method,” J. Acoust. Soc. Am. 113(4), 2077–2083. Yang, A., Lohscheller, J., Berry, D., Becker, S., Eysholt, U., Voigt, D., and Dollinger, M. (2010). “Biomechanical modeling of the three-dimensional J. Acoust. Soc. Am., Vol. 130, No. 1, July 2011 aspects of human vocal fold dynamics,” J. Acoust. Soc. Am. 127(2), 1014–1031. Zemlin, W. R. (1998). Speech and Hearing Science: Anatomy and Physiology, 3rd ed. (Prentice Hall, Englewood Cliffs), pp. 610. Zhang, C., Zhao, W., Frankel, S. H., and Mongeau, L. (2002). “Computational aeroacoustics of phonation, part II: Effects of flow parameters and ventricular folds,” J. Acoust. Soc. Am. 112(5), 2147–2154. Zhang, Z., Neubauer, J., and Berry, D. A. (2007). “Physical mechanisms of phonation onset: a linear stability analysis of aeroelastic continuum model of phonation,” J. Acoust. Soc. Am. 122(4), 2279–2295. Zhao, W., Zhang, C., Frankel, S. H., and Mongeau, L. (2002). “Computational aeroacoustics of phonation, part I: computational methods and sound generation mechanisms,” J. Acoust. Soc. Am. 112(5), 2134–2146. Zheng, X., Bielamowicz, S., Luo, H., and Mittal, R. (2009). “A computational study of the effect of false vocal folds on glottal flow and vocal folds vibration during phonation,” Ann. Biomed. Eng. 37(3), 625–642. Zheng, X. (2009). “Biomechanical modeling of glottal aerodynamics and vocal fold vibration during phonation,” Ph.D. dissertation, George Washington University, DC. Zheng, X., Mittal, R. and Bielamowicz, S. (2011). “A Computational Study of Asymmetric Glottal Jet Deflection During Phonation,” J. Acoust. Soc. Am, 129(4), 2133–2143. Zheng et al.: Three-dimensional laryngeal model 415