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

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 𝑀_∞=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]𝜕𝐖𝜕𝑡+𝜕𝜕𝐅𝑥+𝜕𝜕𝐆𝑦=R1e𝜕𝜕𝐑𝑥+𝜕𝜕𝐒𝑦,(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)𝑒−12𝜌𝑢²+𝑣².(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⁻¹, density 𝜌̂_∞=1.225kgm⁻³, reference length 𝐿_𝑟=0.02m, and dynamic viscosity 𝜂̂_∞=18⋅10⁻⁶Pa⋅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⁻¹ taking into account the tracheal diameter in humans in the range 14.5–17.6 mm [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 𝑔_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 𝐷₁;(2)downstream conditions: 𝑝₂=1/𝜅 and (𝜌,𝜌𝑢,𝜌𝑣) are extrapolated from 𝐷₁;(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𝑘𝐖^𝑛𝑘−1Re𝐑^𝑛𝑘Δ𝑦^𝑘−𝐆^𝑛𝑘−𝑠_2𝑘𝐖^𝑛𝑘−1Re𝐒̃^𝑛𝑘Δ𝑥^𝑘,𝐖^{𝑛+1𝑖,𝑗}=𝜇^{𝑛𝑖,𝑗}𝜇^{𝑛+1𝑖,𝑗}12𝐖^{𝑛𝑖,𝑗}+𝐖^{𝑛+1/2𝑖,𝑗}−Δ𝑡2𝜇^{4𝑛+1𝑖,𝑗}_𝑘=1𝐅^𝑘𝑛+1/2−𝑠_1𝑘𝐖^𝑘𝑛+1/2−1R𝑒𝐑^𝑘𝑛+1/2Δ𝑦^𝑘−𝐆^𝑘𝑛+1/2−𝑠_2𝑘𝐖^𝑘𝑛+1/2−1Re𝐒̃^𝑘𝑛+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 𝐬_𝑘=(𝑠₁,𝑠₂)_𝑘 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Δ𝑡≤CFL||𝑢^max||+𝑐Δ𝑥_min+||𝑣^max||+𝑐Δ𝑦_min+2ReΔ𝑥1^2min+1Δ𝑦^2min⁻¹,(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 𝐷₁ in the narrowest channel cross-section at the middle position of the gap.

(a)

(b)

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⁻¹), Reynolds number Re_∞=4481, atmospheric pressure 𝑝₂=1/𝜅 (𝑝₂̂=102942 𝑃𝑎) at the outlet, and wall oscillation frequency 𝑓=100Hz. The computational domain contained 450×100 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 𝑀_max=0.177 (corresponding to the dimension velocity 𝑢̂_max=60.7ms⁻¹). 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 𝐿2 norm of momentum residuals (ρu)

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 𝐷₁ 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 𝑀_max=0.535 (𝑢̂_max=183.5ms⁻¹) 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.

(a) 𝑡=30ms, 𝑔=1.6mm, 𝑀max=0.159(54.5ms−1)

(b) 𝑡=32.5ms, 𝑔=0.4mm, 𝑀max=0.236(80.9ms−1)

(c) 𝑡=35ms, 𝑔=1.6mm, 𝑀max=0.370(126.9ms−1)

(d) 𝑡=37.5ms, 𝑔=2.8mm, 𝑀max=0.097(33.3ms−1)

(e) 𝑡=40ms, 𝑔=1.6mm, 𝑀max=0.162(55.6ms−1)

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=𝑝−𝑝₂. 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 𝑓₀=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 𝑓₀=100, 𝑓₁=550, 𝑓₂=1150, and 𝑓₃=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 𝑓₁=550Hz (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 (𝐹₁=𝑐̂_∞/(4𝐿⋅𝐿_𝑟)=343ms⁻¹/(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(𝑊_𝑖,𝑗)^𝑛).

(a) The prescribed oscillation of gap width 𝑓=100 Hz.

(b) The acoustic pressure in the gap and DFT spectrum—𝑓0=100Hz.

(c) The acoustic pressure at the outlet and DFT spectrum—𝑓0=100, 𝑓1=550, 𝑓2=1150, 𝑓3=1950Hz.

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ă 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.