Evaluating the accuracy of one-dimensional glottal flow model in predicting voice production: comparison to experiments and three-dimensional flow simulations

A. Vocal fold model

The three-dimensional vibratory dynamics of the vocal folds were simulated using transient response analysis based on mode superposition. This approach significantly reduces computational cost by considering only low-order vibration modes, while still capturing the characteristic vibratory behavior of the vocal folds²⁸.

The vocal folds were modeled as a linear viscoelastic continuum and discretized using the finite element method (FEM). The equation of motion for each node can be expressed as:

where, $\mathbf{M},\mathbf{C}$, and ,$\mathbf{K}$ are mass, damping, and stiffness matrices, respectively, and $\mathit{\delta }$ is the nodal displacement vector. The external force vector $\mathbf{F}$ includes aerodynamic pressure and contact forces between the vocal folds.

In the mode superposition, the nodal displacement $\mathit{\delta }$ is represented as a linear combination of eigenmode ${\mathit{\varphi }}_i$ and modal displacements $q_i$ in the modal space:

Here, $i$ denotes the mode number, and the summation is taken up to the $N\mathrm{th}$ mode. The eigenmodes ${\mathit{\varphi }}_i$ and the corresponding natural angular frequencies ${\omega }_i$ are obtained by solving the characteristic equation derived from the undamped, unforced (in vacuo) system equation:

Solving this eigenvalue problem yields the set of mode shapes ${\mathit{\varphi }}_i$ and their associated frequencies ${\omega }_i$. In general, high-order modes typically lie beyond the frequency range of interest and have negligible influence on voice production. Thus, only a limited number of low-order modes are retained to reduce the number of degrees of freedom in the simulation. Following previous work²⁸, we considered up to the 20th mode in this study, which is the minimal number of modes required to obtain reasonable predictions, resulting in eigenmodes up to 200 Hz being included in the simulation. However, more recent studies²⁹ have reported that considering up to 100 modes may provide a more accurate representation of the detailed vibratory characteristics of the vocal folds (see details in Appendix A).

By expressing the displacement vector as a superposition of eigenmodes, the equation of motion in Eq. (1) can be diagonalized due to the orthogonality of the eigenmode matrix, yielding:

Here, $m_i,c_i$, and $k_i$ are the modal mass, modal damping, and modal stiffness, respectively, and the external force is given by $f_i={\mathit{\varphi }}_i^{T}F$. The angular frequency satisfies the relationship $k_i={\omega }_i^2{m}_i$. Normalizing each eigenmode with respect to the modal mass as ${\tilde{\mathit{\varphi }}}_i=\sqrt{m_i}{\mathit{\varphi }}_i$, the force term as ${\tilde{f}}_i={\tilde{\mathit{\varphi }}}_i^{T}F$, and defining the modal damping ratio as ${\zeta }_i=c_i/\left(2\sqrt{m_i{k}_i}\right)$, Eq. (4) becomes:

By specifying the damping ratio ${\zeta }_i$, the modal displacement $q_i$ can be computed for each mode. The nodal displacement at each time step can then be reconstructed using Eq. (2) by superposing the eigenmodes. It should be noted that this method assumes linear superposition of modal displacements, and therefore does not account for nonlinear modal interactions that may occur during large-amplitude vibrations. However, nonlinear effects induced by external forces, such as subharmonics and transition to chaos, can still be represented²⁹. In this study, we assumed a constant damping ratio of $\zeta =0.06$ for all modes and performed time integration of Eq. (5) using a fourth-order Runge-Kutta method.

The geometry of the vocal folds and vocal tract was simplified as shown in Fig. 1. In this study, the x, y, and z axes corresponded to the inferior–superior, medial–lateral, and anterior–posterior directions, respectively. The vocal tract and subglottal airway were modeled as rectangular ducts. Their cross-sectional area was set to 17 × 16.8 mm², and their lengths were 175 mm and 150 mm, respectively, based on typical adult male anatomy³⁰. An inlet chamber, representing the lungs, was placed upstream of the subglottal airway, and the steady inflow conditions were applied at the entrance.

The vocal fold geometry is based on the widely used M5 model³¹. The geometric parameters r₁ and r₂, as shown in the figure, were set to 6 mm and 0.987 mm, respectively, following the experimental setup³². The width of the vocal fold at the lateral wall was 10.5 mm. The glottal length L_g was set to 17 mm, corresponding to the typical length in adult males³⁰. The vocal fold was placed between the vocal tract and the subglottal airway. The lateral and anterior-posterior surfaces that were attached to the walls were modeled as fully fixed.

The material properties and geometrical parameters used in this study are summarized in Table I. The stiffness and mass matrices were computed based on experimentally measured material properties (Young’s modulus, density, and Poisson’s ratio) of the silicone resin used in previous experiments³³. The eigenvalue analysis was conducted using COMSOL Multiphysics 6.1 (COMSOL Inc.). The vocal fold geometry was discretized using 5650 finite elements. The contact force between the vocal folds was computed using the penalty method³⁴.

B. 3D flow model

In the 3D flow model, the airflow around the vocal folds was predicted by solving the Navier–Stokes equations. To account for both flow and sound generation, we solved the governing equations for a compressible fluid using the finite-difference method. To handle the movement of the vocal fold surface within a structured grid, the immersed boundary method was employed. Specifically, we adopted the volume penalization (VP) method³⁵, in which a penalty term was added as an external force to model the effect of vocal fold surface on the glottal flow.

The governing equations are given by

where $\mathbf{Q}$ represents the conservative variables, $\mathbf{F},\mathbf{G},\mathbf{H}$ are the inviscid flux vectors, and ${\mathbf{F}}_{\nu },{\mathbf{G}}_{\nu },{\mathbf{H}}_{\nu }$ are the viscous flux vectors. The external force $\mathbf{V}$ corresponds to the penalty term in the VP method, and a mask function $\chi$ was used to distinguish between fluid and wall regions. $x_i$ denotes $x_i=(x,y,z)$. In regions where $\chi =1$, the medium was treated as a porous medium. By setting porosity at $\varphi =0.25$, the acoustic reflectivity of the walls was confirmed to be above 99%³⁶.

The mask function near the moving wall was computed as

where d denotes the distance from the wall surface to the computational grid point, and Δx is the grid spacing. By defining the mask function to vary smoothly with distance from the wall, the wall surface can move smoothly across the grid points³⁷.

To accurately compute both the airflow and acoustic pressure fields, a sixth-order accuracy compact scheme was used for spatial differentiation. In addition, a tenth-order accuracy spatial filter was applied to perform implicit large eddy simulation (LES), which models the viscosity of turbulent vortices at the subgrid scale³⁸. For time integration, a third-order accuracy Runge–Kutta method was employed. Details of the methodology are described in the previous study³⁹.

Since the governing equations for airflow were solved using the finite difference method, a structured computational grid was constructed throughout the entire computational domain. Figure 2 shows the computational mesh near the vocal folds. The minimum grid spacing in this region was 0.025 mm, which has been shown to be sufficient for resolving the airflow near the vocal folds in a prior grid convergence study²³. The grid spacing was increased with distance from the vocal folds to reduce computational cost. The total number of grid points used in the fluid simulation was approximately 1.55 × 10⁸. The mesh used in the finite element modeling of the vocal folds was relatively coarse compared to that of the fluid simulation to reduce computational costs. Nevertheless, we confirmed that this mesh resolution had negligible effect on the overall vibratory characteristics (see Appendix B).

Figure 3 illustrates the boundary conditions used in the 3D flow model. In the inlet chamber simulating the lungs, a constant pressure condition was applied, with p_lung = 1000 Pa simulating normal phonation and p_lung = 1800 Pa simulating loud phonation⁴⁰. While the glottal flow region had a relatively fine grid resolution, the vocal tract region had a coarser grid to reduce computational cost. Nevertheless, the maximum grid spacing in the vocal tract was set to 4.86 mm so that sound wave propagation up to approximately 10 kHz can be resolved with at least eight points per wavelength.

A buffer zone was placed outside the acoustic region to prevent reflected waves from reentering the vocal tract. In this zone, the grid spacing was further increased, and a damping force was applied to attenuate the outgoing waves⁴¹. At the outermost boundaries of the domain, non-reflecting and periodic boundary conditions were imposed.

At each time step, the pressure on the vocal fold surface was used to compute the aerodynamic force. Based on this force, the vocal fold displacement and velocity were updated using Eq. (2) and (5). These updated surface displacement and velocity values were then imposed on the flow field. The time integration was carried out by repeating this procedure. The time step for both the flow and vocal fold oscillation solvers was set to 0.25 × 10⁻⁷ s, for both the flow and vocal fold models, ensuring that the Courant number remained below 0.4 based on the speed of sound at 20°C (= 343 m/s).

C. 1D flow model

The 1D flow model was computed using an equivalent electrical circuit model, as described in detail by Ishizaka and Flanagan³, in which components such as resistors, inductors, and capacitors are used to represent the flow resistance, inertia, and compressibility, respectively. The cross-sectional area $A$, circumference $S$, and length $l$ of the subglottal tract and vocal tract were used to calculate the inductance $L=\rho l/2A$, resistance $R=\alpha \left(Sl/A^2\right)\sqrt{\rho \mu \omega /2}$, and capacitance $C=lA/\rho c^2$, where $\rho$ is the air density, $\omega$ is the natural frequency, and $c$ is the speed of sound. The resistance was set in the subglottal tract, and the scaling factor of the resistance $\alpha$ was adjusted to $\alpha =4.00$ under $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$ and $\alpha =2.94$ under $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$, consistent with the previous study²³. Here, $\alpha$ represents the pressure resistance at the sudden contraction from the inlet pressure chamber to the subglottal tract and is tuned to match the amplitude of the subglottal pressure.

The pressure on the vocal fold surface was computed at discrete points from the inferior to the superior point (on the x-axis) using the subglottal pressure p₀, based on Bernoulli’s principle with an added viscous loss term⁷

Here, $U_g$ is the volume flow rate through the glottis, $a_i$ is the cross-sectional area at the $i\mathrm{th}$ point in the medial-lateral $(y-z)$ plane, $\mu$ is the air viscosity, and $\mathrm{\Delta }x$ is the distance between discrete points. $h_{\mathit{eq}}$ represents the height difference between adjacent points and was computed as $h_{eq}=\left(a_i+a_{i+1}\right)/2L_g$. A total of 22 points were placed on the surface of the vocal folds to calculate the pressure distribution.

For simplicity, the flow separation point on the vocal fold surface was assumed to be located at the position of the minimum glottal area. The pressure up to this point was calculated using Eq. (9), while the pressure in the downstream region was assumed to be constant and no pressure recovery was considered.

At the inlet, the lung pressure p_lung was set to 1000 Pa and 1800 Pa as an input voltage in the same way as the 3D flow model. At the outlet of the vocal tract, the transmission line was terminated by a radiation load of a duct with an infinite baffle. At each time step, the flow rate was computed based on the minimum glottal area and the upstream and downstream pressures in the equivalent circuit model. The pressure on the vocal fold surface was then calculated using Equation (9), from which the aerodynamic force acting on each point was obtained. The vocal fold displacement was updated by solving Equation (2) and (5) using this force. This process was repeated at each time step. The time step size was set to 1.0 × 10⁻⁵ s. For the 3D flow model, the simulation was performed on a supercomputer using 100 CPU nodes in parallel, which required approximately 450 hours for each case. In contrast, the 1D flow model can be computed on a single CPU of a laptop in about 8 seconds.

A. Comparison with experimental measurements

We first evaluated to what extent the computational vocal fold model replicated the vocal fold vibration pattern observed in a silicone vocal fold model experiment. We compared the simulation results with the superior view images of the vibrating silicone vocal folds recorded using a high-speed camera³². In this study, the computational model was constructed to have the same geometry and material properties as those in the experiment. For this comparison, the 1D flow model was used, and no vocal tract was included as in the experiment.

Figure 4 shows the superior view images of the vocal folds during one oscillation cycle as observed in both the experiment and the simulation. The fundamental frequency f₀ in the experiment was approximately 133 Hz, while f₀ was about 131 Hz in the simulation. Thus, in both models, one vibration cycle lasted around 7.5 ms. In both sets of images, time zero corresponds to the moment of minimum glottal opening area, and after 3 ms, the vocal folds reached their maximum glottal opening. At this time, the vocal fold vibration amplitude was larger in the middle compared to the anterior and posterior regions. The mean closed quotient (CQ) was 0.18 in the experimental measurements, whereas CQ was 0.11 in the simulation. The glottal closure duration of both the experiment and simulation was slightly shorter than that typically observed in human vocal folds, which is probably due to the relatively small vocal fold thickness⁴² and the use of isotropic materials, which were different from the anisotropic characteristics of actual vocal fold tissue.

Although the large vertical motion of the vocal fold medial edge observed at 8.3 ms in the experiment was not reproduced in the computational model, the shape of the glottal opening exhibited similar trends in all images. This shows that, as in previous studies¹⁸, the eigenmode-based vocal fold model of this study was able to model vocal fold vibrations with sufficient accuracy. The discrepancies between the computation and experiment, such as the shorter glottal closure duration and the absence of large vertical edge motions, might be reduced by refining the grid resolution, as shown in Appendix B. Furthermore, accounting for the effect of geometric nonlinearity associated with large vibration amplitude, as well as nonlinear modal interactions, could also improve the results.

It should be noted that the current 1D flow model was mainly validated for normal, symmetric phonation. In clinical and pathological conditions involving asymmetric or irregular vocal fold vibrations, three-dimensional effects such as complex flow separation, non-uniform tissue deformation, and interactions between the vocal folds may become more significant.

B. 3D flow and pressure fields

In the following, we examine the flow and pressure fields within one vibration cycle. The 3D flow field in the central plane of the vocal folds (z = 0) at a lung pressure of 1800 Pa is shown in Fig. 5. The figure presents the flow filed in different phases, from complete glottal closure at $t/T=0$, through the opening phases at $t/T=0.25\text{and}0.5$, to the closing phase at $t/T=0.75$.

At $t/T=0$, the jet flow from the vocal folds ceased, and turbulent air remained downstream of the glottis. As the vocal folds began to open at $t/T=0.25$, the turbulent airflow was developed and directed straight towards the downstream region of the vocal tract. By the time $t/T=0.75$ in the closing phase, the maximum flow velocity reached $70\mathrm{m}/\mathrm{s}$. In this case, the Reynolds number calculated from the mean flow rate and the glottal length $L_{\mathrm{g}}$ was $Re=U_{\mathrm{g}}/\left(L_{\mathrm{g}}\nu \right)=1848$. The flow structures produced in the 3D model were similar to those observed in previous particle image velocimetry measurements². In that study, a laminar core region extended about 5 mm downstream from the vocal folds, followed by a transitional region of approximately 10 mm, and the jet became fully turbulent around 15 mm downstream. In addition, vortices were formed in the jet shear layer, which were consistent with our observations.

Figure 6 compares the surface pressure distribution obtained from the 3D model at each time step with the 1D model. The horizontal axis of the surface pressure distribution in the 3D flow (Fig. 6 (b)) corresponds to the points P0 to P21 indicated in Fig. 6 (a). Figure 6 (c) compares the surface pressure distribution observed in the 1D model with that obtained at three cross-sections in the 3D model.

At time $t/T=0$, although the vocal folds were almost closed, a small flow rate was observed in both the 1D and 3D models (see Fig. 9 (b) below). As a result, in the 1D model, the pressure remained nearly constant from the subglottal region to the glottal exit, whereas in the 3D model, there was a noticeable difference between the subglottal and supraglottal pressures. The discrepancy at the downstream surface (> P18) indicates that the 1D vocal tract model has difficulty in predicting the downstream pressure when the vocal folds were almost closed.

As the glottis began to open, the surface pressure on the vocal folds gradually decreased towards the downstream direction. Regarding this pressure drop at $t/T=0.25$, the 3D and 1D flow models showed good agreement. Although the 3D model exhibited some differences in pressure along the anterior-posterior direction, these variations were minor (< 200 Pa) compared to the overall pressure drop from the subglottal region (> 1000 Pa).

At $t/T=0.5$, when the glottis was fully open, the subglottal pressure were similar between the 1D and 3D flow models. However, the 1D flow model underestimated the pressure downstream of the glottis due to the overestimated flow rate (see Fig. 9). During the closing phase at $t/T=0.75$, the three-dimensional model showed considerable pressure differences between the center and the anterior or posterior ends of the glottis. Although the 1D flow model fell within the range of these pressures, the 1D flow pressure could not capture such variations, which likely contributed to differences in displacement of the vocal folds.

C. Comparison of 1D and 3D flow models

Figure 7 plots the displacement of the vocal fold center (x = y = z = 0) under $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}\text{and}1800\mathrm{Pa}$. The time is normalized by the vibration period, with t = 0 corresponding to the complete closure at the middle of oscillation cycles. The fundamental frequency of the 3D flow model was 120.8 Hz at $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$ and 125.4 Hz at $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$, while in the 1D flow model, f₀ was 122.6 Hz at $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$ and 127.7 Hz at $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$. In both flow models, the waveform varied slightly from cycle to cycle, resulting in quasi-periodic vibrations. For $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$, the mean difference was 0.17 mm in the x-direction and 0.07 mm in the y-direction, and for $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$, the mean difference was 0.46 mm in the x-direction and 0.10 mm in the y-direction.

In terms of the overall waveform, the displacement of the 1D flow model was slightly smaller than that of the 3D model for $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$, whereas the displacement of the 1D flow model was slightly larger for $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$. Compared to our previous study using a two-mass model²³, the present study with 3D vocal fold vibration exhibited a smaller discrepancy in the medial–lateral (y-axis) displacement (~0.10 mm vs. ~0.42 mm). In the 3D flow model at the higher subglottal pressure, the y-axis displacement reached zero in every cycle, indicating complete glottal closure, whereas in the 1D model, the displacement reached zero only intermittently. The flow separation point predicted by the 3D flow model was located slightly downstream of the minimum glottal area, whereas the 1D flow model assumed separation at the minimum glottal area (as shown in Fig. 6 (c)). In addition, the pressure recovery downstream from the separation point slightly increased the supraglottal pressure values in the 3D flow, which was not included in the 1D flow model. These discrepancies might have contributed to the difference in the vibration amplitude.

Figure 8 displays the subglottal pressure waveforms measured just below the glottis (at x = −10 mm in the 3D model). For both lung pressures of $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$ and $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$, the pressure fluctuation was slightly larger in the 1D flow model. The mean pressure differences were 32 Pa and 160 Pa for $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$ and $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$, respectively. The discrepancy between the 1D and 3D models might be attributed to variations in the glottal opening waveform, as shown in Fig. 7.

At $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$, small ripples appeared in both models, similar to those observed in our previous study²³, due to subglottal and/or vocal tract acoustics. The mean subglottal pressures were 859.5 Pa at $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}\text{and}1527\mathrm{Pa}$ at $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$ in the 1D flow model, while in the 3D flow model, the mean subglottal pressures were 847.2 Pa at $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}\text{and}1428\mathrm{Pa}$ at $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$, respectively. This difference may also have contributed to discrepancies in vocal fold vibration amplitudes.

Figure 9 presents the waveform of the airflow passing through the glottis and the time derivative of the flow rate, dU_g/dt, which represents the intensity of monopole sound sources. According to Lighthill’s acoustic analogy⁴³, the time derivative of the mass flow rate can be regarded as a monopole sound source, and the magnitude of this derivative is often used as an indicator of the source strength¹. To examine the frequency characteristics of the source strength, the spectra of dU_g/dt are also plotted at the bottom.

At $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$, the maximum flow rate was approximately 500 cm³/s, and the 1D flow slightly overestimated the flow rate by approximately 70 cm³/s. At $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$, the maximum flow rate of the 1D flow model was 1330 cm³/s, which was larger than that of the 3D model (970 cm³/s). The mean flow rate difference was 72 cm³/s for $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$ and 185 cm³/s for $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$. The overestimation of the flow rate in the 1D flow model originates from the overestimation of the upstream pressure shown in Fig. 8.

The overall waveforms of dU_g/dt were similar between the 1D and 3D flow models. At $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$, the small ripples observed in the 3D model were partially reproduced by the 1D model, while the 1D flow model overestimated the source strength by 0.5 m³/s² compared to the 3D flow model. In contrast, at $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$, the 1D flow model produced source strengths similar to those of the 3D model, with the maximum difference being approximately 0.21 m³/s². These results suggest that, under lung pressures typical of normal phonation, the 1D flow model can reproduce the results of the 3D flow model with good accuracy.

The spectrum of dU_g/dt was calculated using the discretized Fourier transform with a Hann window. The spectrum was normalized by its maximum value to focus on the comparison in the spectral shapes between the 1D and 3D models. Because the waveforms of dU_g/dt were similar between the 1D and 3D flow models, the overall spectral shapes were also similar. Compared with the 3D flow model, the 1D flow model exhibited slightly higher spectral amplitudes at frequencies from 360 Hz to 840 Hz at $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$, and at 390 Hz at $p_{\text{lung}}=1800\mathrm{P}\mathrm{a}$. In contrast, at $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$, the 1D flow model underestimated the source amplitude above 1600 Hz.

The spectra of pressure at the mouth opening (x = 175 mm) are shown in Fig. 10. The spectral amplitudes were normalized so that the maximum amplitude was set to 0 dB. The first three acoustic resonances (formant frequencies) of the vocal tract were approximately 480 Hz, 1420 Hz, and 2380 Hz, based on the vocal tract length (175 mm) and a speed of sound of 340 m/s. For both $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}1800\mathrm{Pa}$, the 3D and 1D flow models showed spectral peaks at these resonance frequencies. For $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$, the amplitudes at frequencies from 500 to 700 Hz in the 1D flow were approximately 5 dB larger than those in the 3D flow, and these frequencies matched the overestimated peaks observed in the source spectra in Fig. 9. In contrast, at $p_{\text{lung}}=1000\mathrm{P}\mathrm{a}$, the 1D flow model underestimated the source components above 1600 Hz, leading to weaker higher-frequency harmonics at the lips. The relatively stronger white noise band observed in the 3D flow model suggests that the 1D model cannot capture the turbulence-induced noise generated by incomplete glottal closure in the 3D simulation.