Acta Polytechnica

This paper is interested in the mathematical modelling of the voice production process. The main attention is on the possible closure of the glottis, which is included in the model with the concept of a fictitious porous media and using the Hertz impact force The time dependent computational domain is treated with the aid of the Arbitrary Lagrangian-Eulerian method and the fluid motion is described by the incompressible Navier-Stokes equations coupled to structural dynamics. In order to overcome the instability caused by the dominating convection due to high Reynolds numbers, stabilization procedures are applied and numerically analyzed for a simplified problem. The possible distortion of the computational mesh is considered. Numerical results are shown.


Introduction
The voice production mechanism is a complex process consisting of a fluid-structure-acoustic interaction problem, where the coupling between fluid flow, viscoelastic tissue deformation and acoustics is crucial, see [1]. The so-called phonation onset (flutter instability) for a certain airflow rate with a certain prephonatory position leads to the vocal folds oscillation. The important aspect of the phenomena is the glottis closure (glottis is the narrowest part between the vibrating vocal folds).
The considered problem can be mathematically described as a problem of fluid-structure interaction with the involvement of the (periodical) contact problem of the vocal folds. In order to include the interactions of the fluid flow with solid body deformation an the contact problem, a simplified model problem is considered. This model is similar to the simplified two mass model of the vocal folds of [2], see also the aeroelastic model in [3].
In this paper the mathematical model is introduced and the numerical approximation of the problem is described, where the residual based stabilization is used for the incompressible flow model. The simplified lumped vocal fold model is considered, where the contact is treated with the aid of the Hertz impact forces. In the flow model, the contact is considered with the combination of a suitable modification of the inlet boundary condition and the concept of a fictitious porous media flow. Attention is paid to the details of the finite element approximations with the aid of the ALE conservative method. The solution of the system of equations is discussed, with attention on the projection method and on the discrete projection method. Numerical tests are presented.

Mathematical model
In this section a simplified model fluid-structure problem is considered. The domain occupied by fluid Ω t is shown in Figure 1. The fluid flow is coupled with the elastic structure deformation of a simplified two degrees of freedom model of a vocal fold.

Flow problem
First, the air flow is modelled by the system of the Navier-Stokes equations (cf. [4]) written in the ALE conservative form (cf. [5]) where u = (v 1 , v 2 ) is the fluid velocity vector, ρ is the constant fluid density, τ f = (τ f ij ) is the fluid stress tensor given by τ f = −pI + 2µD, the symmetric gradient tensor D = (d ij ) is given by D(u) = 1 2 (∇u + ∇ T u), p denotes the pressure and µ > 0 is the constant fluid viscosity. Further, by A t the Arbitrary Lagrangian-Eulerian mapping is denoted which maps at any time t ∈ [0, T ] the reference configuration Ω ref = Ω 0 onto the current configuration Ω t , by J A the Jacobian of this mapping is denoted, w D denotes the domain velocity, and D A u Dt is the ALE derivative, i.e. the derivative with respect to the reference configuration Ω ref .
For the system (1) the initial and boundary conditions are prescribed. The boundary conditions are prescribed on the boundary ∂Ω f t of the computational domain formed by mutually disjoint parts where Γ I denotes the inlet, Γ O the outlet, Γ S the axis of symmetry and Γ W t denotes either the fixed or the deformable wall. The standard boundary conditions are prescribed at Figure 1. The computational domain Ωt with specification of the boundary parts.
as Γ S is chosen to be located at line y = const (in computations we choose y = 0). Furthermore, at the inlet Γ I and at the outlet part Γ O the boundary conditions are formally prescribed as where n denotes the unit outward normal vector to ∂Ω f t , u I is a prescribed inlet velocity, p ref is a reference pressure value (p ref = 0 in what follows), ε > 0 is a penalization parameter and α − denotes the negative part of a real number α. Here, the boundary condition (3c) weakly imposes the Dirichlet boundary condition u = u I with the aid of a penalization parameter ε.

Structure model
The motion of the vocal fold is modelled as a motion of a rigid body governed by the displacements w 1 (t) and w 2 (t) of the two masses m 1 and m 2 , respectively (see Fig. 2). The equation of motion (see [3] for details) reads Mẅ where M is the mass matrix of the system, K is the diagonal stiffness matrix of the system characterized by spring constants k 1 , k 2 , and B = ε 1 M + ε 2 K, is the matrix of the proportional structural damping, ε 1 , ε 2 are the constants of the proportional damping. The mass matrix is given by

Coupling conditions
The aerodynamical forces F 1 , F 2 are computed with the aid of the aerodynamical lift force L(t) and the aerodynamical torsional moment M (t) acting on the surface of the structure Γ W t . The aerodynamical lift force and the aerodynamical torsional moment are evaluated with the aid of the mean (kinematic) pressure p and the mean flow velocity u = (u 1 , u 2 ) as the integrals over the surface of the airfoil where l denotes the depth of the profile section, and the vector r ort has components r ort being the position of the structure elastic axis.
The displacement of the structure surface Γ W t is determined in terms of w 1 , w 2 , and it is used as the boundary condition for the displacement of any point of the computational domain Ω ref 0 onto the domain Ω t . Moreover, the domain velocity at Γ W t is determined using the w 1 (t) and w 2 (t) obtained from solution of the ordinary differential equations (4).

Treatment of the contact problem
In order to close of the fluid computational domain, the idea of a fictitious porous media flow is employed. It means that during the closing phase of the vocal folds, the part of the structure domain (denoted by Ω p t ) is assumed to be a domain of porous media flow, see Fig. 3. The flow in the domain Ω p t is then assumed to be governed by the modified equations where the coefficient α corresponds to the artificial porosity of the medium, see [6]. The coefficient in the porous media flow is usually chosen as α = µ P where µ is the dynamical viscosity of the fluid and P is an artificial permeability coefficient, see [6]. In practical computations equation (7) is solved in the whole computational domain Ω t with α being zero at all points outside of the porous media domain Ω P t and being a positive constant (in this paper α was chosen as 10 7 ) in the interior of Ω P t . The porous media domain is determined by the following procedure illustrated in Figure 3. First, the actual gap g(t) at time t is computed as distance of Γ W t from the axis of symmetry Γ S and compared to a prescribed minimal gap threshold g min . If g(t) ≥ g min , then the glottis is open and no porous media is used (i.e. Ω P t is empty and α ≡ 0). For g(t) < g min the y-displacement of the points of the surface of Γ W t is modified in order to not violate the condition g(t) ≥ g min , i.e. these points are vertically displaced by such a shift −∆w, which makes the actual gap equal to g(t) = g min . For the structural model the Hertz impact forces are used as specified in [3]. In this case, the right hand side F of Eq. (4) is modified by the addition of the Hertz impact force F H . The Hertz impact force is then distributed to the both components of F depending on the position x max of the impact. The model of Herz impact forces is given as where δ stands for the penetration of the vocal fold through the contact plane, b H is a damping factor (here set to zero), r is the radius of the osculating circle (i.e. inverse of the curvature), E is Young's modulus of elasticity of the vocal fold and µ H is its Poisson's ratio, for values see [3].

Mesh deformation
First, due to the motion of the structure, the transformation of the computational domain (and mesh as well) needs to be constructes at every time instant t k .
Here it is represented by the construction of the ALE mapping A t at every time instant t = t k . The ALE mapping A t is sought in the form of displacements d = d(ξ, t) of points ξ of the original configuration, i.e. A t (ξ) = ξ + d(ξ, t). The boundary condition for d is then known for every point of the boundary ∂Ω: at ∂Ω \ Γ W t the displacement d = 0 and at Γ W t it is determined by the already known (or estimated) position of the structure surface Γ W t in terms of the displacements w 1 and w 2 . In order to determine d in Ω a boundary value problem is solved, see e.g. [7].

Numerical approximation
In this section let us consider the computational domain Ω t and the ALE mapping A t to be already known and A t to be sufficiently smooth at every time instant t ∈ I = (0, T ). A similar assumption is made about the domain velocity w D .

Weak formulation
In order to introduce the function spaces for velocity and pressure, we start with defining the spaces for velocity and pressure at the reference configuration. We shall seek the velocity-pressure pair U = (u, p) at Further, the weak formulation is derived using the standard approach, i.e. we multiply equations (1) Next, these (reference) test functions are defined on the current configuration Ω t with the aid of the ALE mapping A t , i.e. the test function z = z(x, t) can be defined using a reference test function We shall denote the space of such test functions by X , i.e.
The space X is for any t subspace of V t , and the test functions from X are time independent when transformed on the reference domain. This means that their ALE derivative equals zero, i.e. for z ∈ X we have In what follows by the symbol (·, ·) M the dot product in L 2 (M) or L 2 (M) is denoted. The weak form is then derived in the standard form: we take test function V = (z, q) ∈ X × Q t multiply the first equation (1) by test functions z and the second by q, sum together, integrate over Ω t and use Green's theorem for viscous terms and the pressure gradient. Thus we arrive to the weak form: Find U = (u, p) ∈ W t such that for any t ∈ (0, T ) u satisfy the boundary condition (2a) and holds for any V = (z, q) ∈ X × Q t . The forms in Eq. (11) are defined for any U = (u, p) ∈ W t , U = (u, p) ∈ W t and V = (z, q) ∈ X × Q t as follows: the form d in Eq. (11) is defined by the boundary forms b, L Γ are given as and the convective term is given by the skewsymmetric trilinear form c (here we abbreviate w The time derivative terms arises from the identity which follows from (see [8]) In what follows an equivalent integral form of equation (13) shall be used in the form for any volume V t whose motion is fully determined by the ALE mapping, i.e. V t = A t (V ref 0 ).

Time discretization
For the purpose of time discretization an equidistant partition t j = j∆t of the time interval I with the constant time step ∆t > 0 is considered and we denote the approximations of velocity and pressure by u j ≈ u(·, t j ) and p j ≈ p(·, t j ) for j = 0, 1, . . . , u j ∈ V tj and p j ∈. Similarly, by J A j and w j D the Jacobian of the ALE mapping A t at time instant t j and the domain velocity at time instant t j is denoted. In what follows we shall describe the discretization at a fixed time instant t n+1 . The time derivative on the right hand side of equation (12) is approximated at t = t n+1 formally by the second order backward difference formula or more precisely the time derivative term from (11) is approximated by where the test function z ∈ X is a time dependent function defined by (9)

Finite element approximations
The spaces V tn+1 and X are approximated using the finite element subspaces V h and X h constructed over an admissible triangulation T ∆ of the domain Ω, respectively. Similarly, the pressure space Q tn+1 is approximated by its finite element subspace Q h constructed again over T ∆ . Here, the Taylor-Hood finite elements are used, i.e. the spaces of continuous piecewise quadratic velocities and continuous piecewise linear pressures are used, i.e. the velocities are sought in the space of the test functions is given by and the trial/test pressure space is given as vol.

Special Issue/ Numerical solution of FSI problems with contacts
The finite element approximations of u and p are then sought in the finite element spaces V h = X h × Q h constructed over an admissible triangulation τ h of the computational domain Ω f t : Find an approximate solution U h = (u, p) ∈ W h such that Eq. (19) holds for any test function V h = (z, q) ∈ X h ×Q h . Furthermore, this formulation is stabilized using the SUPG/PSPG stabilization terms together with the div-div stabilization terms given as where the modified test function is given by Φ(U ; V ) = ((u − w D ) · ∇)z + ∇q and δ K , τ K are suitably chosen stabilization parameters, see e.g. The stabilized discrete formulation then reads: Find holds for any test function

Linearization
In order to solve the nonlinear problem (21), the sequence of the linearized problem is solved until it converges to a sufficient precision. We start with an approximation U 0 and for k = 0, 1, . . . solve the linearized problems: Find U =: U k+1 such that holds for any V = (z, q) ∈ X h × Q h . Using the finite element base functions the discretization leads to the system of linear equations in the form of a saddle point problem, i.e. in the form where the matrix A = M + D + C + A B + A S consists of the mass matrix M, the diffusion matrix D, the convection matrix C, the boundary terms matrix A B and the stabilization matrix A S part arising from the forms M , d, c, B and S, respectively. Matrix B * is the discrete gradient operator. It corresponds to the form b * including the stabilization terms, B T corresponds to the discrete divergence operator realized by the form b plus stabilization. The matrix Q follows fully from the stabilization S, i.e. in particular for the case with no stabilization the matrix Q = 0. The solution of such a problem is difficult, see e.g. [9]. This is caused by the presence of the continuity equation, which can be treated at the continuous level by the approach based on the projection method. On the other hand, this can be understood as an approximation of the solution of the system (23), or more precisely as preconditioned iterative solution of this system.

Projection method
The projection method is based on the Helmholtz-Hodge decomposition of any vector field v, i.e. v = v div + ∇Ψ, where v div is the divergence free field, see [10]. In this section, for the sake of simplicity, we shall restrict ourselves to the case of the fixed computational domain Ω and to the case of the first order discretization in time. This means that here we shall discuss the solution of the Navier-Stokes equations in the form where the time derivative term is replaced by the backward difference formula Here, the problem (24) can be sought using the segregated approach, see [11]. This can be formally written as I. Solve momentum equations forũ II.Projectũ on the space of (discrete) divergence free functions: Solve pressure equation III.Update the velocity field by adding the pressure gradient: Let us mention that the steps I. and II. need to be equipped with suitable boundary conditions. Due to splitting the coupled equations, the boundary conditions as presented in (3-2) require a modification. Particularly the boundary conditions (2) contain pressure, which is not available in the step I, but can be used from the previous time step. The pressure equation needs to be equipped by the Neumann type boundary condition, where the compatibility condition needs to be satisfied for the existence of a solution, which is unique to a constant. Nevertheless, due to splitting these two steps, the splitting error arises, for overview see e.g. [12]. In the considered problem, the main difficulty is in the realization of the non-standard boundary conditions (2,3). To this end we shall consider another possibility based on the preconditioned solution of the discrete problem (23).

Discrete projection methods and preconditioning
Let us consider the system of linear equations written in the form where the matrix M is given as This is the well studied case, where the matrix A is (possibly) a non-symmetric positive definite and the matrix B has full rank due to the Babuška-Brezzi inf-sup condition satisfied. The matrix M can be factorized as where S is the pressure Schur complement given as The inverse matrix to the matrix M given (28) is able to compute the inverse matrix to the matrix A and the inverse to the Schur complement matrix S. In order to solve the problem (27), the three sub-steps can be followed similarly as in Section 4.1. In this process, the solution of the system with matrix A and the matrix S needs to be realized. Here, matrix A is a mass matrix perturbed with convectiondiffusion. The stabilization terms and the solution of the system with the matrix A can be realized effectively. However the solution of the system with the Schur complement matrix S must be realized iteratively. This is why the above approach is used only as a preconditioner for GMRES method, where the matrix S is replaced by a suitable approximation, see e.g. [13]. In the considered discretization, the matrix of the system is given by where the matrix A is (possibly) a non-symmetric positive definite, the matrix B has full rank and matrix D is positive semidefinite. In this case the matrix M can be factorized as where S 1 is the pressure Schur complement given as

Numerical Results
First, the Oseen problem is considered in the computational domain Ω = (0, 1) 2 . The problem is equipped with the Dirichlet boundary condition u = b prescribed at ∂Ω. Here, α = 0 is used, b = (sin(πx), −πy cos(πx)) and the right hand side f is chosen in such a way that u = (sin(πx), −πy cos(πx)) is solution of the Oseen problem. The computations were performed for different values of the viscosity coefficient nu. First, the convergence of the Galerkin finite element approximations u G h to the exact solution u = b is investigated, p(x, y) = sin(πx) cos(πy) for ν = 0.05 (here, relatively high viscosity was chosen in order to obtain stable Galerkin approximations even on coarser meshes). For an approximation of the flow problem, the Taylor Hood finite elements were used. The errors in the H 1 norm are shown in Table 1, where by h max the length of the maximal edge in the triangulation is denoted. These results are compared to the results of the stabilized formulation of the same problem, which shows that the used residual based stabilization does not pollute the solution, see Table 2. The convergence orders in both cases agree well with the theoretical estimate. For the stabilized method such convergence rates are well preserved for the values ν = 10 −3 , . . . , 10 −6 with a slow down observed only for coarse grid configuration.     Let us emphasize, that as the convection b equals the exact solution u, the Dirichlet problem for Navier-Stokes equations can be formulated with the same analytical solution. This problem was solved as an Navier-Stokes problem with a known analytical solution. One can really notice that b = u, and thus the given u is the analytical solution of the Navier-Stokes problem with the known Dirichlet boundary condition. Similarly as above, the numerical approximation was compared for both the stabilized and non-stabilized method, with the results confirming the expected theoretical order of convergence in H 1 -norm, see Table 3 for the stabilized method.
Further, the method for the solution of FSI problem with contact was realized and tested on a benchmark problem from [14]. The results for the inflow velocity U ∞ = 0.65 m/s are shown in terms of displacements w 1 and w 2 in dependence on time in Figures 5-6. One can easily see that the system becomes aeroelastically unstable and the so called phonation onset phenomena arises. With further continuation, the vibrating vocal folds start to be influenced by their mutual contact, see the vibrations resembling a limit cycle oscillations in Figure 6. Similar behaviour is observed also for the inflow velocity U ∞ = 0.7 m/s with a much faster appearance of the contact. Figure 9 then shows the details of the flow in terms of flow velocity in the glottis region.

Conclusion
In this paper the mathematical description of a simplified fluid-structure interaction problem arising in the biomechanics of voice production was presented. Special attention was paid to considering the contact of the vibrating vocal folds, which was treated with the aid of the Hertz impact forces, a suitable choice of inlet boundary conditions and the use of a fictitious porous media flow description. Attention was also given to the numerical discretization of the problem, to the solution of the linearized problem and to the realization of the gap closing. The numerical tests were shown to verify the theoretical error estimates of the applied method and the numerical results of a benchmark test problem were presented.