Numerical approximations of flow induced vibrations of vocal folds

The paper focus on mathematical modelling of incompressible fluid flow interacting with vibrations of an elastic vocal fold. The flow in moving domain is modelled by the incompressible Navier-Stokes equations written in the Arbitrary Lagrangian-Eulerian (ALE) form. The channel geometry is an approximation of the human glottal region. The flow model is coupled with a simplified structure model. The problem is mathematically described and the resulting fluid-structure interaction problem is discretized by a stabilized finite element method. A strong coupling algorithm is applied for solution of the coupled fluid-structure problem. The choice of boundary conditions is discussed, particularly the choice of different artificial inlet/outlet boundary conditions is described in details. The numerical results are shown.


Introduction
The mathematical modelling of fluid-structure interactions is a complex difficult task. Nowadays it has become to be important also in biomechanics of human voice, see [1]. Naturally in medicine it is important to understand the basic principles of voice production, see e.g. [2]. However, the voice production is a complex mechanism which includes complex airflow, vocal fold vibrations, acoustic resonances, etc. The voice sound is created in the glottis -the narrow part of the channel between the two vocal folds. The vocal folds vibrations leads to periodical (almost or complete) closure of the channel. The collision of the vocal folds creates the sound, but mathematical and numerical modelling of such a complex phenomenon encounters many difficulties as e.g. the correct choice of boundary conditions at the artificial inlet/outlet parts of boundary or the treatment of the mesh deformation.
One possibility how to study the complex problem of voice creation is the use of experimental methods. The experiments partly focus on measurement of the externally driven models, e.g. [3], [4], [5]), or on self-oscillating models, e.g. [6], [7], [8]. However, there are many difficulties encountered in both in vivo and in vitro measurements.
The mathematical models on the other hand can be used to provide important information about the voice production mechanisms. In this case several simplifications are usually made for the purpose of the computational models, for example see the two-mass model of the vocal folds, see [9] or several works using the same concept as [10]. These simple mass-spring models coupled with a quasi-1D airflow proved to be very useful, see e.g. [10]), and they can be used to estimate the frequena Corresponding author: petr.svacek@fs.cvut.cz cies or vocal folds loading by impact stress, see [11]. The main advantage of these models is the ability to perform nearly real-time simulations. However, these mathematical models usually do not address the complex problem in all details and uses several simplifications.
Particularly, the simplified flow models do not describe the glottal airflow details. That is why recently the numerical solution of the 2D Navier-Stokes equations has been used. For instance see [12], where the finite volume approximations of the Navier-Stokes equations were coupled with the two-mass dynamic model. See also [13], where the 2D finite element approximations of the Navier-Stokes equations coupled with a two degrees of model with the vibrating vocal fold described by a mechanically equivalent two degrees of freedom system. The fluid-structure interaction using the coupled finite element (FE) models of the vocal folds and airflow was published in [6], [14]. Nevertheless also in these high-fidelity models several simplifications are used, e.g. no contact problem is considered.
In this paper we shall focus on addressing the proper inlet/outlet boundary conditions applicable in the case of closing of the channel due to the vocal folds vibrations. First, the mathematical model addressing the fluid structure interactions is described. Further, the numerical approximation of the described mathematical model by the finite element method is given. The numerical results are presented and the choice of boundary conditions and its implementation is discussed. To this end the Arbitrary Lagrangian Eulerian method is applied in order to treat the motion of the computational domain. Further the motion of the vocal folds is described as linear 2D elastic continuum on the reference domain, i.e. Lagrangian concept is applied. In order to simplify the model a two degrees of freedom structure is used. The interface conditions are presented.

Arbitrary Lagrangian Eulerian method
First, the time dependence of the computational fluid flow domain Ω f t needs to be addressed. Here, the socalled arbitrary Lagrangian-Eulerian (ALE) is used, see e.g. [15]. In order to apply this method a regular smooth one-to-one ALE mapping A of the reference configuration Ω f 0 onto the current configuration Ω f t , i.e it maps the point X from the reference configuration X ∈ Ω f 0 on the point x of the current configuration x ∈ Ω f t at any time instant t. The point x with a fixed reference X moves in time x = A(X, t) with the domain velocity defined by Further, the time derivative with respect to the reference frame is introduced as ALE derivative. The ALE derivate is then given as where the functionf is the transformation of the function f onto the reference frame, i.e.f (X, t) = f (A(X, t), t) for any X ∈ Ω f 0 . The ALE derivative then satisfy (see [15])

Incompressible Navier-Stokes equations
The flow of an incompressible viscous fluid in the domain Ω f t is described by the system of the incompressible Navier-Stokes equations (see, e.g., [16]) written in the ALE form Here is the fluid velocity vector, ρ is the constant fluid density, and τ f = (τ f ij ) is the fluid stress tensor given by where p is the pressure, μ > 0 is the constant fluid viscosity, δ ij denotes the Kronecker symbol and D is the symmetric part of the velocity gradient given by For the system (4) the initial condition v(x, 0) = v 0 (x) for x ∈ Ω f 0 is specified. Further the boundary conditions needs to be prescribed on the boundary ∂Ω f t of the computational domain. The boundary ∂Ω f t is assumed to be formed by mutually disjoint parts where Γ I denotes the inlet part of the boundary, Γ O is the outlet part of the boundary and Γ W t denotes the (fixed and moving) wall. The following boundary conditions are prescribed where p ref denotes a reference mean pressure value at the outlet (we set p ref = 0) and n is the unit outward normal vector to ∂Ω t and v I is the inlet velocity. For the analysis of the boundary condition 7c) see e.g. [17], [18]. Particularly this type of boundary condition treats well the cases of the (backward) inflow of the fluid at the outlet part of the boundary Γ O . For this reason this boundary condition can be also used at the inlet boundary. This means that the condition 7, a) can be replaced by the condition The boundary condition (8)  In this paper both approaches were tested.

Structure model
Elastic structure. The deformation of vocal folds is described as deformation of elastic structure. Let us consider a bounded domain Ω s ⊂ R 2 representing the vocal fold. The displacement of the structure is denoted by u, which is a function u = u(X, t) defined at a time instant t ∈ (0, T ) and a point X = (X 1 , X 2 ) ∈ Ω s . The deformation of the elastic body is modelled by the dynamic elasticity. The equations of motion of an elastic body have the form where f = (f 1 , f 2 ) is the density of the volume force, ρ s denotes the density of the structure and the stress tensor τ s is given by the generalized Hooke law for isotropic material in the form τ s = λ s e kk (u) I + 2μ s e(u)x where λ s and μ s are the so-called Lamé coefficients. The small strain tensor e = (e ij ) then is given by where ∇ denotes the gradient operator with respect to variables X 1 , X 2 . Similarly as for fluid, the initial and boundary conditions are prescribed for system (9, i.e. the initial conditions read u(·, 0) = 0 and ∂u ∂t (·, 0) = 0. The boundary ∂Ω s is then decomposed in two disjoint parts, the moving interface  Reduced order model. The vibrations of the vocal folds usually exhibit two fundamantal frequencies. Such vibrations can be described with an simplified aeroelastic two degrees of freedom model. The motion of Γ W t is governed by the displacements θ 1 (t) and θ 2 (t) (upward positive) of the two masses m 1 and m 2 , respectively (see Fig. 3). The displacement vector θ = (θ 1 , θ 2 ) T is obtained by the solution of the following equations (see [10] for details) where M is the mass matrix, K is the diagonal stifness matrix with spring constants c 1 , c 2 on its diagonal and B is the matrix of the proportional structural damping. The mass matrix is given by where m 1 , m 2 , m 3 are masses. The components of F = (F 1 , F 2 ) T are the aerodynamical forces (downward positive). The proportional damping matrix is chosen as The aerodynamic forces F 1 F 2 in Eq. (13) depends on flow velocity v and pressure p.
where L is the length of the vocal fold, h is the depth of the structure (i.e. third dimension in z axis), value p f (x, t) denotes the value of the pressure at time instant t on the surface of the vocal fold positioned at x and n 2 (x, t) denotes the second component of the unit outward normal to the surface of the vocal fold at position x and time t.

Numerical approximation
The presented fluid-structure problem is first decoupled. The motion equations are time integrated with the aid of Newmark method. The fluid flow is approximated by the stabilized finite element method. The coupled fluidstructure model is solved with the aid of strongly coupled scheme.
The main attention is paid to the approximation of the fluid flow model. However, the application of the finite element method for incompressible Navier-Stokes equations requires to carefully treat several issues. The first one is the dominating advection in the momentum equations. The second is the compatibility of the velocity/pressure spaces. Moreover, in order to obtain physically admissible solution it is usually necessary to apply suitable mesh refinement (e.g. anisotropically refined mesh, see [19]) combined with a stabilization technique, see [20], [21], [22]. The applied stabilization scheme needs to be modified in the case of the application on moving domains (see [23]).
First the problem is time discretized, i.e the time interval is divided 0 = t 0 < t 1 < · · · < T, t k = kΔt with a constant time step Δt > 0. The velocity v and the pressure p at time t n are then approximated by v n and p n , respectively. Furthermore, the ALE velocity w D (t n+1 ) is approximated by w n+1 D computed at a point x = A(X, t n+1 ) by (15) The ALE derivative of the velocity v at time t n+1 is approximated using the second-order backward difference formula, i.e.
where at a given time instant t n+1 we denotev The time discretized problem is then formulated weakly with the aid of the velocity space W = H 1 (Ω), the pressure space Q = L 2 (Ω) and the space of test functions Now, the weak formulation of the time discretized problem (4) is obtained with the multiplication by the test functions ϕ ∈ X,q ∈ Q, sumation, integration over Ω f t and transformation with the aid of Green's theorem. The weak solution U = (v, p) then satisfy equation where the form a(U, V ) is defined by and the convective form c(v * , v, ϕ) reads is obtained from the convective term (v * · ∇v, ϕ) Ω by integration by parts and using the boundary conditions.
In order to apply the Galerkin finite element method (FEM) to the discretization of the problem 18, we approximate the spaces W , X, Q by their finite dimensional subspaces W h , X h , Q h , i.e.
The couple (X h , Q h ) of the finite element spaces is chosen based on a regular triangulation of the computational domain such that it satisfy the Babuška-Brezzi (BB) infsup condition (see, e.g., [24], [25]). The spaces W h , X h and Q h are formed by the well-known Taylor-Hood P 2 /P 1 conforming finite elements used for the velocity/pressure approximation. This means that p h is a linear function and v h is a quadratic vector-valued function on each element K ∈ T h . The fully stabilized scheme is used (see, e.g., [26], [21]). The stabilized discrete problem reads: The stabilization terms are defined using the modified test function defined on each element K as . The stabilization terms then reads as and the div-div stabilization reads Here, δ K ≥ 0 and τ K ≥ 0 are suitable chosen parameters. The choice of the parameters δ K and τ K is carried out according to [26] and [22]. The parameter δ K is defined on the basis of the local transport velocity, local element size and local viscozity. In the case of the Taylor-Hood finite elements and a bounded convection v ≈ 1 the following choice of parameters appears suitable: where τ * > 0 and δ * > 0 are fixed constants.

Numerical results
For the computations the fluid density ρ = 1.2 kg m −3 , and fluid kinematic viscosity ν = 1.58 × 10 −5 m 2 /s were chosen. The distance of the masses from the center was l = L/2, the lengths of sub-and supra-glottal regions were L 0 = 2L and L 2 = 2L, respectively. The height of the channel was 2H 0 , where H 0 = 1 2 g 0 + max x∈ 0,L a(x) and g 0 is the initial gap g(0) = g 0 . The initial gap g 0 was chosen 0.4 mm and 0.6 mm, respectively.
Using the given shape and dimension of the vocal fold and using the structural density ρ h the total mass m, the moment of inertia I and the excentricity e is computed. The system is then replaced by the equivalent three mass system where m 1,2 = 1 2l 2 (I + m e 2 ± m e l) and m = m 1 + m 2 + m 3 . The structural parameters were chosen according [27], see Tab. 1. The presented results shows the computed case with the prescribed inlet velocity. The prescribed inlet pressure did not lead to an instability in the considered range of velocities/pressures.
The solution of the aeroelastic system was performed for two cases represented by Tab. 1. The solution of the flow model (4), structure model (13) and interface conditions was numerically approximated in a simplified geometry shown in Figs. 1, 2. The aeroelastic response for the is shown in Fig. 4 for Model F and in Fig. 5 for Model M in time domain in terms of displacements θ 1 (t) and θ 2 (t). The vibrations of the structure die in time after a time period due to both structural and strong aerodynamic damping for lower flow velocities. With increasing flow velocity the flutter type of instability can be observed. The critical velocity is in a good agreement with results presented in [27]. The physical meaning of the instability is so-called phonation onset which is an important voice production characteristic in humans. The flow velocity patterns are shown in Fig. 6, where the typical pressure distribution and flow velocity magnitude pattern can be observed.

Conclusion
In this paper a fluid structure interaction problem of flow induced vibrations of vocal folds was described. The numerical approximations based on the finite method were presented. The choice of the inlet and the outlet boundary condition was discussed. On the considered model problem only the prescribed inlet velocity formulations were successful. The prescribed pressure difference lead to an additional non-physical damping of the resulting aeroelastic system. This is particularly caused by the fact, that the prescribed inlet pressure boundary condition is imposed only weakly, and during the computation the pressure is not constant in time, but depends on the vocal fold vibrations.