Abstract
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 , Reynolds number 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] 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 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 , density , reference length , and dynamic viscosity . Inflow temperature [] depends on relation for speed of sound where 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 . The nondimensional dynamic viscosity in the dissipative terms is a function of temperature in the form . The heat transfer coefficient is expressed as , where 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โ taking into account the tracheal diameter in humans in the range 14.5โ17.6โ [2].
The bounded computational domain , 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 and maximum , not closing the channel completely.
The boundary conditions are considered in the following formulation:(1)upstream conditions: , , , and is extrapolated from domain ;(2)downstream conditions: and are extrapolated from ;(3)flow on the wall: where is velocity vector of the wall and 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 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 at initial time to a domain at [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] where is the time step, is the volume of cell , and are the steps of the grid in directions and , and vector 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.
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 [8]. Artificial dissipation is used to stabilize computation and also due to velocity gradients in the narrowest width of the channel, where . The vector of conservative variables is computed at a new time level : .
The stability condition of the scheme (on the regular orthogonal grid) limits the time step where denotes the local speed of sound, and are the maximum velocities in the domain, and 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 to capture the boundary layer effects. Figure 3 shows the detail of the grid in domain in the narrowest channel cross-section at the middle position of the gap.
(a)
(b)
5. Numerical Results
The numerical results were obtained (using a specifically developed program) for the following input data: Mach number (), Reynolds number , atmospheric pressure (โ) at the outlet, and wall oscillation frequency . The computational domain contained cells in .
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 and the convergence to the steady-state solution computed using the norm of momentum residuals (). The maximum Mach number computed in the flow field Figure 4(a) was (corresponding to the dimension velocity ). 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.
(a) The velocity flow field mapped by isolines of Mach number
(b) Convergence to the steady-state solution using norm of momentum residuals (ฯu) |
The numerical simulation of the airflow computed in domain 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โโ effect is apparent in the flow field pattern. The absolute maximum of Mach number () in the flow field during fourth cycle was achieved at time โms (โ, opening phase) behind the narrowest channel cross-section. The flow becomes practically periodic after the first period of oscillations.
(a) , , |
(b) , , |
(c) , , |
(d) , , |
(e) ms, , |
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 . 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 โ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 , , , and โHz can be identified in the spectrum envelope of the pressure at the channel outlet in Figure 6(c)-right. The first acoustic resonance (see spectrum peak at cc 500โHz) corresponds to the first eigenfrequency of a simple tube of the length of the complete channel closed at the inlet and open at the outlet (). 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 ).
(a) The prescribed oscillation of gap width โHz. |
(b) The acoustic pressure in the gap and DFT spectrumโ. |
(c) The acoustic pressure at the outlet and DFT spectrumโ, , , . |
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 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โ 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.
Acknowledgment
This paper was partially supported by Research Plans MSM 6840770010, GAฤR P101/11/0207 and 201/08/0012.