Large-Eddy Simulation of Internal Flow through Human Vocal Folds

The phonatory process occurs when air is expelled from the lungs through the glottis and the pressure drop causes flow-induced oscillations of the vocal folds. The flow fields created in phonation are highly unsteady and the coherent vortex structures are also generated. For accuracy it is essential to compute on humanlike computational domain and appropriate mathematical model. The work deals with numerical simulation of air flow within the space between plicae vocales and plicae vestibulares. In addition to the dynamic width of the rima glottidis, where the sound is generated, there are lateral ventriculus laryngis and sacculus laryngis included in the computational domain as well. The paper presents the results from OpenFOAM which are obtained with a large-eddy simulation using second-order finite volume discretization of incompressible Navier-Stokes equations. Large-eddy simulations with different subgrid scale models are executed on structured mesh. In these cases are used only the subgrid scale models which model turbulence via turbulent viscosity and Boussinesq approximation in subglottal and supraglottal area in larynx.


Introduction
Large-eddy simulation was proposed in as early as 1963 by Smagorinsky [7]. In LES the large scale motions of turbulent flow are computed directly and small scale motions (sub-grid scale = SGS) are modelled. LES has significant reduction in computational cost compared to DNS and is more accurate than the RANS approach. The large eddies contain most of the turbulent energy and these eddies are responsible for most of the momentum transfer and turbulent mixing. The small eddies incline to be isotropic and homogenous (and therefore the SGS motions should be easier to model than model all scales within a domain). Nowadays, LES is promising numerical tool for simulating realistic turbulent flows. For instance inside the complicated measurable place in human larynx.

Mathematical formulation
The Navier-Stokes equations are derived from the conservation laws for mass, momentum and energy. In LES a low-pass spatial filter is applied to the instantaneous conservation equations to formulate equations for large scale motions (explicit filtering). If the finite volume method (FVM) is used for solving the instantaneous governing equations numerically, then there is no need to apply a filter to the equations explicitly (it is called implicit filtering, equivalent to convolution with a top-hat filter). The filtered equations in a Newtonian incompressible flow can be written as whereū i is filtered velocity, ρ is density, µ is molecular viscosity, S i j is resolved scale strain rate tensor. [10] In the context of LES τ i j is called the subgrid scale (SGS) Reynolds stress. It plays a role in LES similar to the Reynolds stress in RANS, but physics that it models is different. The SGS energy is usually a much smaller part of the total energy at the turbulent velocity field than the RANS turbulent energy. Thus a model accuracy may be less crucial in LES in comparison with RANS. [1] To model SGS stress tensor we can use eddy-viscosity assumption (Boussinesq's hypothesis) as follows where µ t is SGS eddy viscosity. If µ t is substitute into the equation (2) then becomes (6) and the remaining problem now is how to determine µ t . The Smagorinsky model contains where C S is the Smagorinsky constant (0.18). [4] The value gives reasonable results for isotropic turbulence and for flows near a solid wall it should be 0.1. [10] Nevertheless, this model has shortcomings as for example being too dissipative. Improvement of this model was suggested by Germano [3]. In OpenFOAM (OF) the µ t is determined as follows where C k and C E are model constants (C k = 0.094, C E = 1.048). Symbol ∆ represents a filter width, c is deltaCoeff (in OF=1) and V C is the volume of the cell.

Geometry and mesh
The geometry of the vocal folds was specified by Scherer [6]. For LES computations, a structural mesh from 400k volume cells was created and it included 5 cells in the zaxis (5x 0.2 mm). The one of the important parameter for testing dynamic mesh (where the vocal folds have twodegrees of freedom) is maximal non-orthogonality factor. The analysis of this parametes make sense in the time, when the vocal folds are in the nearest position to each other, i.e. t = 0.003 s and the maximal non-orthogonality factor is 55.14). Mesh spacings in wall units are commonly used to indicate LES adequacy. Especially the theoretical limits for LES and DNS are 10 ≤ ∆x + ≤ 20, ∆y + wall < 1, 5 ≤ ∆z + ≤ 10, where the x is streamwise direction, the y direction is wall normal and the z is spanwise and homogeneous. [2] In this study the y + for the critical time when the mesh is maximally deformed is: y + avg = 1.17, y + max = 17.29 and in different time is usually about y + avg = 0.31, y + max = 3.29.

Numerical solution
In this LES study, central differencing scheme (CDS) is used for spatial discretization of the diffusive term. The CDS can be obtained by a Taylor series expansion where the terms involve derivatives of the second order (higher are neglected). The second order CDS is non-dissipative and conservative which is essential for LES. Upwindbased schemes are not used in LES, because of the pro-duction of high numerical dissipation. For spatial discretization of convective term, a total variation diminishing scheme (TVD) is used. A numerical method is said to be TVD, if the total variation (TV) is not increased in time, see (17), where φ is a variable, i is an index of the node. A monotone scheme is TVD and a TVD scheme is monotonicity preserving. It does not create any new local extrema within solution domain. In our case we have used high resolution scheme, see (18), where the term called the anti-diffusive flux creates 2nd order of accuracy (which it decreases numerical diffusion, but also it leads to unphysical oscillations). To preserve the good properties, i.e. the stability from the 1st order scheme and the accuracy from the 2nd order scheme is to multiple the flux limiter ψ to the anti-diffusive term, see (19), where indexes are Ddownwind (left), C-center, U-upwind (right), f -face (between C and U) and r f is taken as the ratio of two consecutive gradients. [5] In this study was used the flux limiter of the MUSCL TVD scheme as the discretization scheme of the convective term, see (20) and CD TVD scheme with ψ(r f ) = 1 as discretization scheme of the diffusive term.
For the temporal discretization a second-order backward implicit Euler was used. Since the time steps are small in LES, it is not essential to use different temporal scheme than second-order.

Results
After one period which is composed from one convergent and one divergent motion of the vocal folds, this state is captured on Fig. 2. The velocity field shows the airflow from subglottal area to supraglottal area (vocal folds are in divergent position). An air jet in glottal area is formed on the surface of the right vocal fold margin. The velocity of an airflow is much lower (circa 50 % lower) than in the similar study [9] with no turbulent model, because the effective viscosity (ν e f f = ν + ν S GS ) is computed higher in the high shear regions. The velocity is also slightly lower with 2D Smagorinsky and 2D One-Equation compared with [9]. Detailed studies deal with vocal folds are also [11] and [8]. On the Fig. 3 the high shear regions are shown: 1) The shear layer within rima glottidis, 2) The upper parts of the sacculus laryngis walls, 3) The peripherals of large eddies in supraglottal area. In the ideal case the influence of the mesh the mesh should not be seen.  Fig. 4 are shown isosurfaces of the Q-criterion which detect the regions where coherent structures can exist. If the laplacian of the pressure is positive then the condition about generating the vortex tubes is satisfied. In the rima glottidis we can recognize the so-called hairpin horseshoe vortex which is generally formed near the wall as a result of the Tollmien-Schlichting instability. The leg of the vortex is in a buffer layer in the longitudinal direction, the neck and the head of the vortex are also seen. A regular shape of this kind of vortex is usually an exception, because the transition of the turbulence regime is stormy dynamic process in which the shear flow created from the main flow interact with the neighbouring structures. This kind of vortices are significant for generating and self-sustaining the turbulent structures.
The phenomenon in the supraglottal section is the socalled merging of vortices where these vortices have the same orientation of rotation and create cirles. The spatial effect of this mechanism is related to pseudoperiodic disturbances and the energy within the vocal fold domain is transfered from small eddies via merging to create one larger eddy (reverse energy cascade). The interactions of vortex tubes are also captured and it looks like "8". Two vortices are connected in the one point and subsequently disconnected vice versa than how they were connected before (the so-called bridging).

Conclusion
In this study, 3D LES case was simulated with Smagorinsky SGS model on structured mesh in OpenFOAM 4.1. The weaknesses of this study are following. The Smagorinsky model overpredicts the dissipation rate and for further work is highly recommend to use one-equation SGS model with van Driest damping function and more requirement because of the additional transport equation for the SGS kinetic energy. Another point for correct 3D modeling is a guarantee fine mesh refinement for wall-resolved LES in all three directions. Last but not least important point is initialization of turbulence at inlet. Current inflow boundary condition generation methods in LES are generally two basic categories: the precursor method (an addition simulation) and synthesis method (random fluctuations at the inlet). This is the main issue of the further work, because the vocal folds domain is computed from the pressure difference (not by prescribed velocity).