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

Authors: Petra Pořízková; Karel Kozel; Jaromír Horáček
Journal: Journal of Applied Mathematics
 (2012)

---

## 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 numberM∞=0.012, Reynolds number Re∞=4481, and the wall motion frequency 100 Hz.


---

## Body

## 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](2.1)∂W∂t+∂F∂x+∂G∂y=1Re(∂R∂x+∂S∂y),
where W=[ρ,ρu,ρv,e]T is the vector of conservative variables, F and G are the vectors of inviscid fluxes, and R and S are the vectors of viscous fluxes. The static pressure p is expressed by the state equation in the form(2.2)p=(κ-1)[e-12ρ(u2+v2)].
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 ĉ∞=343ms-1, density ρ̂∞=1.225kgm-3, reference length L̂r=0.02m, and dynamic viscosity η̂∞=18·10-6Pa·s. Inflow temperature T̂∞ [K] depends on relation for speed of sound ĉ∞2=κR̂T̂∞ where κ=1.4 is the ratio of the specific heats. Inflow pressure satisfies equation of state for ideal gas p̂∞=ρ̂∞R̂T̂∞.General Reynolds number in (2.1) is computed from reference variables Re=ρ̂∞ĉ∞L̂r/η̂∞. The nondimensional dynamic viscosity in the dissipative terms is a function of temperature in the form η=(T/T∞)3/4. The heat transfer coefficient is expressed as k=ηκ/[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.2ms-1 taking into account the tracheal diameter in humans in the range 14.5–17.6 mm [2].The bounded computational domainD1, 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 g(t) 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 gmin⁡=0.4mm and maximum gmax⁡=2.8mm, not closing the channel completely.Figure 1
Computational domainD1.L=8(160mm),H=0.8(16mm), and g=0.08(1.6mm)—middle position.The boundary conditions are considered in the following formulation:(1)
upstream conditions:u∞=û∞/ĉ∞=M∞, v∞=0, ρ∞=1, and p∞ is extrapolated from domain D1;(2)
downstream conditions:p2=1/κ and (ρ,ρu,ρv) are extrapolated from D1;(3)
flow on the wall:(u,v)=(uwall,vwall) where (uwall,vwall) is velocity vector of the wall and ∂T/∂n⃗=0 where T=κp/ρ is the temperature.The general Reynolds number in (2.1) multiplied with nondimensional value M∞H represents kinematic viscosity scale, and for computation of the real problem, inlet Reynolds number Re∞=ρ̂∞ĉ∞M∞HL̂r/η̂∞ 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 subdomainsDi,j (i.e., quadrilateral cells). The system of (2.1) is integrated over the subdomains Di,j 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 Dt=0 at initial time t=0 to a domain Dt at t>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]
(4.1)Wi,jn+1/2=μi,jnμi,jn+1Wi,jn-Δtμi,jn+1∑k=14[(F̃kn-s1kWkn-1ReR̃kn)Δyk-(G̃kn-s2kWkn-1ReS̃kn)Δxk],W¯i,jn+1=μi,jnμi,jn+112(Wi,jn+Wi,jn+1/2)-Δt2μi,jn+1∑k=14[(F̃kn+1/2-s1kWkn+1/2-1ReR̃kn+1/2)Δyk-(G̃kn+1/2-s2kWkn+1/2-1ReS̃kn+1/2)Δxk],
where Δt=tn+1-tn is the time step, μi,j=∫∫Di,jdxdy is the volume of cell Di,j, Δx and Δy are the steps of the grid in directions x and y, and vector sk=(s1,s2)k represents the speed of edge k (see Figure 2). The physical fluxes F,G,R,S on the edge k of the cell Di,j are replaced by numerical fluxes (marked with tilde) F̃,G̃,R̃,S̃ as approximations of the physical fluxes.Figure 2
Finite volumeDi,j and dual volume Vk′.The approximations of the convective termssWk and the numerical viscous fluxes R̃k,S̃k on the edge k are central. The higher partial derivatives of velocity and temperature in R̃k,S̃k are approximated using dual volumes Vk′ (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 dissipationAD(Wi,j)n [8]. Artificial dissipation is used to stabilize computation and also due to velocity gradients in the narrowest width of the channel, where M∞≈0.5. The vector of conservative variables is computed at a new time level tn+1: Wi,jn+1=W¯i,jn+1+AD(Wi,j)n.The stability condition of the scheme (on the regular orthogonal grid) limits the time step(4.2)Δt≤CFL(|umax⁡|+cΔxmin⁡+|vmax⁡|+cΔymin⁡+2Re(1Δxmin⁡2+1Δymin⁡2))-1,
where c denotes the local speed of sound, umax⁡ and vmax⁡ 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 iny-direction is Δymin⁡≈1/Re∞ to capture the boundary layer effects. Figure 3 shows the detail of the grid in domain D1 in the narrowest channel cross-section at the middle position of the gap.Grid of quadrilateral cells in the narrowest part of domainD1 at the middle position of the gap width g=0.08(1.6mm).Detail: Δymin⁡≐5·10-4(0.01mm).
(a)(b)


## 5. Numerical Results

The numerical results were obtained (using a specifically developed program) for the following input data: Mach numberM∞=0.012 (û∞=4.116ms-1), Reynolds number Re∞=4481, atmospheric pressure p2=1/κ (p̂2=102942 Pa) at the outlet, and wall oscillation frequency f̂=100Hz. The computational domain contained 450×100 cells in D1.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.Figure4 shows the initial condition for unsteady computation of the flow field in domain D1 and the convergence to the steady-state solution computed using the L2 norm of momentum residuals (ρu). The maximum Mach number computed in the flow field Figure 4(a) was Mmax⁡=0.177 (corresponding to the dimension velocity û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.The initial condition computed inD1—M∞=0.012,Re∞=4481,p2=1/κ,450×100 cells, and Mmax⁡=0.177 (ûmax⁡=60.7ms-1).
(a)
The velocity flow field mapped by isolines of Mach number(b)
Convergence to the steady-state solution usingL2 norm of momentum residuals (ρu)The numerical simulation of the airflow computed in domainD1 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 “Coandă’’ effect is apparent in the flow field pattern. The absolute maximum of Mach number Mmax⁡=0.535 (ûmax⁡=183.5ms-1) in the flow field during fourth cycle was achieved at time t=34.2 ms (g=1.002 mm, opening phase) behind the narrowest channel cross-section. The flow becomes practically periodic after the first period of oscillations.The unsteady numerical solution of the airflow inD1—f̂=100Hz,M∞=0.012,Re∞=4481,p2=1/κ,and  450×100 cells. Data computed during the fourth oscillation cycle. Results are mapped by isolines of velocity ratio and by streamlines.
(a)
t=30ms, g=1.6mm, Mmax⁡=0.159(54.5ms-1)(b)
t=32.5ms, g=0.4mm, Mmax⁡=0.236(80.9ms-1)(c)
t=35ms, g=1.6mm, Mmax⁡=0.370(126.9ms-1)(d)
t=37.5ms, g=2.8mm, Mmax⁡=0.097(33.3ms-1)(e)
t=40ms, g=1.6mm, Mmax⁡=0.162(55.6ms-1)Figure6 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 pac=p-p2. 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 f̂0=100 Hz and the generated higher harmonics as a consequence of the throttling and nearly closing the glottal gap g(t). 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 f̂0=100, f̂1=550, f̂2=1150, and f̂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 f̂1=550Hz (see spectrum peak at cc 500 Hz) corresponds to the first eigenfrequency F̂1 of a simple tube of the length of the complete channel closed at the inlet and open at the outlet (F̂1=ĉ∞/(4L·L̂r)=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(Wi,j)n).Three vibration periods of the gap width oscillations and the acoustic pressure signalspac computed in the gap and at the outlet on the axis of the channel.
(a)
The prescribed oscillation of gap widthf̂=100 Hz.(b)
The acoustic pressure in the gap and DFT spectrum—f̂0=100Hz.(c)
The acoustic pressure at the outlet and DFT spectrum—f̂0=100, f̂1=550, f̂2=1150, f̂3=1950Hz.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 “Coandă 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 “Coandă” 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.



---

*Source: 545120-2012-05-13.xml*