Table of Contents Author Guidelines Submit a Manuscript
Journal of Applied Mathematics
VolumeΒ 2012, Article IDΒ 545120, 9 pages
Review Article

Numerical Simulation of Unsteady Compressible Flow in Convergent Channel: Pressure Spectral Analysis

1Department of Technical Mathematics, Faculty of Mechanical Engineering, Czech Technical University in Prague, Karlovo NΓ‘m. 13, 121 35 Praha 2, Czech Republic
2Institute of Thermomechanics AS CR, DolejΕ‘kova 5, 18200 Prague 8, Czech Republic

Received 16 January 2012; Accepted 8 March 2012

Academic Editor: Fu-YunΒ Zhao

Copyright Β© 2012 Petra PoΕ™Γ­zkovΓ‘ et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


This study deals with the numerical solution of a 2D unsteady flow of a compressible viscous fluid in a channel for low inlet airflow velocity. The unsteadiness of the flow is caused by a prescribed periodic motion of a part of the channel wall with large amplitudes, nearly closing the channel during oscillations. The flow is described by the system of Navier-Stokes equations for laminar flows. The numerical solution is implemented using the finite volume method (FVM) and the predictor-corrector Mac-Cormack scheme with Jameson artificial viscosity using a grid of quadrilateral cells. Due to the motion of the grid, the basic system of conservation laws is considered in the arbitrary Lagrangian-Eulerian (ALE) form. The numerical results of unsteady flows in the channel are presented for inlet Mach number π‘€βˆž=0.012, Reynolds number Re∞=4481, and the wall motion frequency 100 Hz.

1. Introduction

A current challenging question is a mathematical and physical description of the mechanism for transforming the airflow energy in human vocal tract (convergent channel) into the acoustic energy representing the voice source in humans. The voice source signal travels from the glottis to the mouth, exciting the acoustic supraglottal spaces, and becomes modified by acoustic resonance properties of the vocal tract [1]. The airflow coming from the lungs causes self-oscillations of the vocal folds, and the glottis completely closes in normal phonation regimes, generating acoustic pressure fluctuations. In this study, the movement of the boundary channel is known, harmonically opening and nearly closing in the narrowest cross-section of the channel, making the investigation of the airflow field in the glottal region possible.

Acoustic wave propagation in the vocal tract is usually modeled from incompressible flow models separately using linear acoustic perturbation theory, the wave equation for the potential flow [2], or the Lighthill approach on sound generated aerodynamically [3].

Goal of this work is numerical simulation of compressible viscous flow in 2D convergent channel which involves attributes of real flow causing acoustic perturbations as is β€œCoandΔƒ phenomenon’’ (the tendency of a fluid jet to be attracted to a nearby surface), vortex convection and diffusion, jet flapping, and so forth along with lower call on computer time, due to later extension in 3D channel flow. Particular attention is paid to the analysis of acoustic pressure signal from the channel.

2. Mathematical Model

To describe the unsteady laminar flow of a compressible viscous fluid in a channel, the 2D system of Navier-Stokes equations was considered as a mathematical model. The system of Navier-Stokes equations is expressed in nondimensional conservative form [4]πœ•π–+πœ•π‘‘πœ•π…+πœ•π‘₯πœ•π†=1πœ•π‘¦ξ‚΅Reπœ•π‘+πœ•π‘₯πœ•π’ξ‚Άπœ•π‘¦,(2.1) where 𝐖=[𝜌,πœŒπ‘’,πœŒπ‘£,𝑒]𝑇 is the vector of conservative variables, 𝐅 and 𝐆 are the vectors of inviscid fluxes, and 𝐑 and 𝐒 are the vectors of viscous fluxes. The static pressure 𝑝 is expressed by the state equation in the form1𝑝=(πœ…βˆ’1)π‘’βˆ’2πœŒξ€·π‘’2+𝑣2ξ€Έξ‚„.(2.2) The transformation to the nondimensional form uses inflow parameters (marked with the infinity subscript) as reference variables (dimensional variables are marked with the hat): the speed of sound Μ‚π‘βˆž=343msβˆ’1, density Μ‚πœŒβˆž=1.225kgmβˆ’3, reference length ξπΏπ‘Ÿ=0.02m, and dynamic viscosity Μ‚πœ‚βˆž=18β‹…10βˆ’6Paβ‹…s. Inflow temperature ξπ‘‡βˆž [𝐾] depends on relation for speed of sound ̂𝑐2βˆžξπ‘…ξπ‘‡=πœ…βˆž where πœ…=1.4 is the ratio of the specific heats. Inflow pressure satisfies equation of state for ideal gas Μ‚π‘βˆž=Μ‚πœŒβˆžξπ‘…ξπ‘‡βˆž.

General Reynolds number in (2.1) is computed from reference variables Re=Μ‚πœŒβˆžΜ‚π‘βˆžξπΏπ‘Ÿ/Μ‚πœ‚βˆž. The nondimensional dynamic viscosity in the dissipative terms is a function of temperature in the form πœ‚=(𝑇/π‘‡βˆž)3/4. The heat transfer coefficient is expressed as π‘˜=πœ‚πœ…/[Pr(πœ…βˆ’1)], where Pr=0.7 is the Prandtl number.

3. Computational Domain and Boundary Conditions

For phonation of vowels, the frequencies of the vocal folds oscillations are in the region from cc 82 Hz for bass up to cc 1170 Hz for soprano in singing voice, and the airflow velocity in the trachea is approximately in the range of 0.3–5.2 msβˆ’1 taking into account the tracheal diameter in humans in the range 14.5–17.6 mm [2].

The bounded computational domain 𝐷1, used for the numerical solution of flow field in the channel, is shown in Figure 1. The domain is a symmetric channel, the shape of which is inspired by the shape [5] of the trachea (inlet part of the channel), vocal folds, false vocal folds and supraglottal spaces (outlet part). The upper and the lower boundaries are the channel walls. A part of the wall changes its shape between the points A and B according to a given harmonic function of time and axial coordinate (see, e.g., [6]). The gap 𝑔(𝑑) is the narrowest part of the channel (in point C). The gap width was oscillating with frequency 100 Hz (typical for normal male voice) between the minimum 𝑔min=0.4mm and maximum 𝑔max=2.8mm, not closing the channel completely.

Figure 1: Computational domain 𝐷1. 𝐿=8(160mm), 𝐻=0.8(16mm), and 𝑔=0.08(1.6mm)β€”middle position.

The boundary conditions are considered in the following formulation:(1)upstream conditions: π‘’βˆž=Μ‚π‘’βˆž/Μ‚π‘βˆž=π‘€βˆž, π‘£βˆž=0, 𝜌∞=1, and π‘βˆž is extrapolated from domain 𝐷1;(2)downstream conditions: 𝑝2=1/πœ… and (𝜌,πœŒπ‘’,πœŒπ‘£) are extrapolated from 𝐷1;(3)flow on the wall: (𝑒,𝑣)=(𝑒wall,𝑣wall) where (𝑒wall,𝑣wall) is velocity vector of the wall and πœ•π‘‡/πœ•βƒ—π‘›=0 where 𝑇=πœ…π‘/𝜌 is the temperature.

The general Reynolds number in (2.1) multiplied with nondimensional value π‘€βˆžπ» represents kinematic viscosity scale, and for computation of the real problem, inlet Reynolds number Re∞=Μ‚πœŒβˆžΜ‚π‘βˆžπ‘€βˆžπ»ξπΏπ‘Ÿ/Μ‚πœ‚βˆž is used.

4. Numerical Solution

The numerical solution uses FVM in conservative cell-centered form on the grid of quadrilateral cells, see, for example, [4].

The bounded domain is divided into mutually disjoint subdomains 𝐷𝑖,𝑗 (i.e., quadrilateral cells). The system of (2.1) is integrated over the subdomains 𝐷𝑖,𝑗 using the Green formula and the mean value theorem. In the time-changing domain, the integral form of FVM is derived using the ALE formulation. The ALE method defines homeomorphic mapping of the reference domain 𝐷𝑑=0 at initial time 𝑑=0 to a domain 𝐷𝑑 at 𝑑>0 [7].

The explicit predictor-corrector MacCormack (MC) scheme in the domain with a moving grid of quadrilateral cells is used. The scheme is 2nd order accurate in time and space using orthogonal regular grid [4] 𝐖𝑛+1/2𝑖,𝑗=πœ‡π‘›π‘–,π‘—πœ‡π‘›+1𝑖,𝑗𝐖𝑛𝑖,π‘—βˆ’Ξ”π‘‘πœ‡4𝑛+1𝑖,π‘—ξ“π‘˜=1ξ‚π…ξ‚ƒξ‚€π‘›π‘˜βˆ’π‘ 1π‘˜π–π‘›π‘˜βˆ’1𝐑Reπ‘›π‘˜ξ‚Ξ”π‘¦π‘˜βˆ’ξ‚€ξ‚π†π‘›π‘˜βˆ’π‘ 2π‘˜π–π‘›π‘˜βˆ’1̃𝐒Reπ‘›π‘˜ξ‚Ξ”π‘₯π‘˜ξ‚„,𝐖𝑛+1𝑖,𝑗=πœ‡π‘›π‘–,π‘—πœ‡π‘›+1𝑖,𝑗12𝐖𝑛𝑖,𝑗+𝐖𝑛+1/2𝑖,π‘—ξ‚βˆ’Ξ”π‘‘2πœ‡4𝑛+1𝑖,π‘—ξ“π‘˜=1ξ‚π…ξ‚ƒξ‚€π‘˜π‘›+1/2βˆ’π‘ 1π‘˜π–π‘˜π‘›+1/2βˆ’1𝐑Rπ‘’π‘˜π‘›+1/2ξ‚Ξ”π‘¦π‘˜βˆ’ξ‚€ξ‚π†π‘˜π‘›+1/2βˆ’π‘ 2π‘˜π–π‘˜π‘›+1/2βˆ’1̃𝐒Reπ‘˜π‘›+1/2Δπ‘₯π‘˜ξ‚„,(4.1) where Δ𝑑=𝑑𝑛+1βˆ’π‘‘π‘› is the time step, πœ‡π‘–,𝑗=βˆ«βˆ«π·π‘–,𝑗𝑑π‘₯𝑑𝑦 is the volume of cell 𝐷𝑖,𝑗, Ξ”π‘₯ and Δ𝑦 are the steps of the grid in directions π‘₯ and 𝑦, and vector π¬π‘˜=(𝑠1,𝑠2)π‘˜ represents the speed of edge π‘˜ (see Figure 2). The physical fluxes 𝐅,𝐆,𝐑,𝐒 on the edge π‘˜ of the cell 𝐷𝑖,𝑗 are replaced by numerical fluxes (marked with tilde) ̃𝐒𝐅,𝐆,𝐑, as approximations of the physical fluxes.

Figure 2: Finite volume 𝐷𝑖,𝑗 and dual volume π‘‰ξ…žπ‘˜.

The approximations of the convective terms π¬π–π‘˜ and the numerical viscous fluxes ξ‚π‘π‘˜,Μƒπ’π‘˜ on the edge π‘˜ are central. The higher partial derivatives of velocity and temperature in ξ‚π‘π‘˜,Μƒπ’π‘˜ are approximated using dual volumes π‘‰ξ…žπ‘˜ (see [4]) shown in Figure 2. The inviscid numerical fluxes are approximated by the physical fluxes from the cell on the left side of the current edge in the predictor step and from the cell on the right side of the current edge in the corrector step.

The last term used in the MC scheme is the Jameson artificial dissipation AD(π‘Šπ‘–,𝑗)𝑛 [8]. Artificial dissipation is used to stabilize computation and also due to velocity gradients in the narrowest width of the channel, where π‘€βˆžβ‰ˆ0.5. The vector of conservative variables is computed at a new time level 𝑑𝑛+1: 𝐖𝑛+1𝑖,𝑗=𝐖𝑛+1𝑖,𝑗+AD(π‘Šπ‘–,𝑗)𝑛.

The stability condition of the scheme (on the regular orthogonal grid) limits the time step||𝑒Δ𝑑≀CFLmax||+𝑐Δπ‘₯min+||𝑣max||+𝑐Δ𝑦min+21ReΞ”π‘₯2min+1Δ𝑦2minξƒͺξƒͺβˆ’1,(4.2) where 𝑐 denotes the local speed of sound, 𝑒max and 𝑣max are the maximum velocities in the domain, and CFL<1 for nonlinear equations [9]. Time discretization of the scheme satisfies discrete geometric conservation law (DGCL), see [10].

The grid used in the channel has successive refinement cells near the wall. The minimum cell size in 𝑦-direction is Δ𝑦minβˆšβ‰ˆ1/Re∞ to capture the boundary layer effects. Figure 3 shows the detail of the grid in domain 𝐷1 in the narrowest channel cross-section at the middle position of the gap.

Figure 3: Grid of quadrilateral cells in the narrowest part of domain 𝐷1 at the middle position of the gap width 𝑔=0.08(1.6mm). Detail: Δ𝑦min≐5β‹…10βˆ’4(0.01mm).

5. Numerical Results

The numerical results were obtained (using a specifically developed program) for the following input data: Mach number π‘€βˆž=0.012 (Μ‚π‘’βˆž=4.116msβˆ’1), Reynolds number Re∞=4481, atmospheric pressure 𝑝2=1/πœ… (̂𝑝2=102942β€‰π‘ƒπ‘Ž) at the outlet, and wall oscillation frequency 𝑓=100Hz. The computational domain contained 450Γ—100 cells in 𝐷1.

The computation has been carried out in two stages. First, a steady numerical solution is obtained, when the channel between points A and B has a rigid wall fixed in the middle position of the gap width. Then this solution is used as the initial condition for the unsteady simulation.

Figure 4 shows the initial condition for unsteady computation of the flow field in domain 𝐷1 and the convergence to the steady-state solution computed using the 𝐿2 norm of momentum residuals (πœŒπ‘’). The maximum Mach number computed in the flow field Figure 4(a) was 𝑀max=0.177 (corresponding to the dimension velocity ̂𝑒max=60.7msβˆ’1). The picture displays nonsymmetric flow developed behind the narrowest channel cross-section. The graph in Figure 4(b) indicates the nonstationary solution of initial condition which is caused probably by eddies separated in the unmovable glottal orifice and floating away.

Figure 4: The initial condition computed in 𝐷1---π‘€βˆž=0.012, Re∞=4481, 𝑝2=1/πœ…, 450Γ—100 cells, and 𝑀max=0.177 (̂𝑒max=60.7msβˆ’1).

The numerical simulation of the airflow computed in domain 𝐷1 over the fourth cycle of the wall oscillation is presented in Figure 5 showing the unsteady flow field in five time instants during one vibration period. Large eddies are developing in supraglottal spaces and a β€œCoandΜ†a’’ effect is apparent in the flow field pattern. The absolute maximum of Mach number 𝑀max=0.535 (̂𝑒max=183.5msβˆ’1) in the flow field during fourth cycle was achieved at time 𝑑=34.2 ms (𝑔=1.002 mm, opening phase) behind the narrowest channel cross-section. The flow becomes practically periodic after the first period of oscillations.

Figure 5: The unsteady numerical solution of the airflow in 𝐷1β€” 𝑓=100Hz, π‘€βˆž=0.012, Re∞=4481, 𝑝2=1/πœ…, and 450Γ—100 cells. Data computed during the fourth oscillation cycle. Results are mapped by isolines of velocity ratio and by streamlines.

Figure 6 shows three vibration periods of the gap width oscillation (a) and the acoustic pressures signals computed in the gap (b) and at the outlet (c) on the axis of the channel. The acoustic pressure was calculated by subtracting the average values of the pressure signals 𝑝ac=π‘βˆ’π‘2. The acoustic pressure time dependent data are transformed to frequency-dependent data (acoustic pressure spectrum) using discrete fourier transformation (DFT) of the signals. The spectrum of the pressure in the glottis, Figure 6(b)-right, shows the dominant fundamental frequency of vocal folds model oscillations 𝑓0=100 Hz and the generated higher harmonics as a consequence of the throttling and nearly closing the glottal gap 𝑔(𝑑). Two different regimes are apparent in the acoustic pressure signal at the channel outlet during one vibration period of the glottis in Figure 6(c)-left. Relatively smooth signal containing low frequencies is dominant in the time interval corresponding to the phase of maximum glottal opening, and a very noisy signal containing high frequencies is associated with the phase of minimum glottal opening. Four acoustic resonances of the channel cavity at about 𝑓0=100, 𝑓1=550, 𝑓2=1150, and 𝑓3=1950 Hz can be identified in the spectrum envelope of the pressure at the channel outlet in Figure 6(c)-right. The first acoustic resonance 𝑓1=550Hz (see spectrum peak at cc 500 Hz) corresponds to the first eigenfrequency 𝐹1 of a simple tube of the length of the complete channel closed at the inlet and open at the outlet (𝐹1=Μ‚π‘βˆžξπΏ/(4πΏβ‹…π‘Ÿ)=343msβˆ’1/(4β‹…0.16m)=536Hz). The acoustic resonances are more damped with increasing frequency, which can be caused by the fluid viscosity as well as by a numerical viscosity implemented in the numerical method (constants magnitude in AD(π‘Šπ‘–,𝑗)𝑛).

Figure 6: Three vibration periods of the gap width oscillations and the acoustic pressure signals 𝑝ac computed in the gap and at the outlet on the axis of the channel.

Remark 5.1. We used several tests on fine and coarse grids in the computational domain. Also the domain has been prolonged in upstream and also in downstream part. Achieved results were approximately the same on fine grid.

Remark 5.2. The mathematical model (2.1) of laminar flow used in this case is debatable. For the first approximation, we supposed unformed turbulent flow at the inlet part of the channel.

Remark 5.3. The validation of computations for this case is not complete because of experiments absence. Semivalidation of the computations is comparison with particle image velocimetry method (PIV) experiment, but we can compare only qualitative behavior of the flow. Full validation of the code for subsonic and transonic flow through a turbine cascade computed in periodic domain is showed, for example, in [10].

6. Discussion and Conclusions

The numerical solution in the channel showed large vortex structures developed in the supraglottal space moving slowly downstream and decaying gradually. It was possible to detect a β€œCoandΜ†a phenomenon” in the computed flow field patterns. A similar generation of large-scale vortices, vortex convection and diffusion, jet flapping, and general flow patterns were experimentally obtained in physical models of the vocal folds by using PIV method in [5, 11, 12].

The results show that some numerical results of viscous flow in a symmetric channel using a symmetric grid and scheme can be nonsymmetrical, depending on the geometry and the Reynolds number. This effect was observed also for laminar transonic flow computation [13]. The assumption of the axisymmetry solution for the axisymmetry channels (see [6]) excludes modeling the β€œCoandΜ†a” effect and large vortex structures of the size comparable with the cross-section of the channel.

The analysis of the computed pressure revealed basic acoustic characteristics of the channel. This is promising result for future studies for a direct modeling of human voice generated by the flow in vibrating glottis taken into account a real shape of the human vocal tract for phonation and throttling the glottal gap width up to zero.


This paper was partially supported by Research Plans MSM 6840770010, GAČR P101/11/0207 and 201/08/0012.


  1. I. R. Titze, Principles of Voice Production, National Centre for Voice and Speech, Iowa City, IA, Iowa, USA, 2000.
  2. I. R. Titze, The Myoelastic Aerodynamic Theory of Phonation, National Centre for Voice and Speech, Iowa City, IA, Iowa, USA, 2006.
  3. S. Zorner, M. Kalteenbacher, W. Mattheus, and C. Brucker, β€œHuman phonation analysis by 3d aero-acoustic computation,” in Proceedings of the International Conference on Acoustic NAG/DAGA 2009, pp. 1730–1732, Rotterdam, The Netherlands, 2009.
  4. J. FΓΌrst, M. Janda, and K. Kozel, β€œFinite volume solution of 2D and 3D Euler and Navier-Stokes equations,” in Mathematical Fluid Mechanics, pp. 173–194, BirkhΓ€user, Basel, Switzerland, 2001. View at Google Scholar Β· View at Zentralblatt MATH
  5. J. Horacek, P. Sidlof, V. Uruba, J. Vesely, V. Radolf, and V. Bula, β€œPIV measurement of flow-patterns in human vocal tract model,” in Proceedings of the the International Conference on Acoustic NAG/DAGA 2009, pp. 1737–1740, Rotterdam, The Netherlands, 2009.
  6. P. PunčochΓ‘Ε™ovΓ‘-PoΕ™Γ­zkovΓ‘, J. FΓΌrst, J. HorÑček, and K. Kozel, β€œNumerical solutions of unsteady flows with low inlet Mach numbers,” Mathematics and Computers in Simulation, vol. 80, no. 8, pp. 1795–1805, 2010. View at Publisher Β· View at Google Scholar Β· View at Zentralblatt MATH
  7. R. Honzatko, J. Horacek, and K. Kozel, Solution of Inviscid Incompressible Flow Over a Vibrating Profile, vol. 3 of COE Lecture Notes, Kyushu university, 2006, Edited by M. Benes, M. Kimura and T. Nataki.
  8. A. Jameson, W. Schmidt, and E. Turkel, β€œNumerical solution of the Euler equations by the finite volume methods using Runge-Kutta time-stepping schemes,” AIAA, pp. 81–125, 1981. View at Google Scholar
  9. R. Dvorak and K. Kozel, Mathematical Modeling in Aerodynamics, Prague, Czech Republic, 1996, (In Czech).
  10. P. Porizkova, Numerical solution of compressible flows using finite volume method, Ph.D. thesis, Czech Technical University in Prague, Faculty of Mechanical Engineering, Prague, Czech Republic, 2009.
  11. P. Sidlof, Fluid-structure interaction in human vocal folds, Ph.D. thesis, Charles University in Prague, Faculty of Mathematics and Physics, Prague, Czech Republic, 2007.
  12. J. Neubauer, Z. Zhang, R. Miraghaie, and D. Berry, β€œCoherent structures of the near field flow in self-oscillating physical model of the vocal folds,” Journal of the Acoustical Society of America, vol. 121, no. 2, pp. 1102–1118, 2007. View at Google Scholar
  13. J. Fort, T. Hulek, K. Kozel, and M. Vavrincova, Remark on Numerical Simulation of 2D Unsteady Transonic Flows, vol. 414 of Lecture Notes in Physics, Springer, Berlin, Germany, 1993, Edited by M. Napolitano and F. Sabetta. View at Publisher Β· View at Google Scholar