Numerical solution of fluid-structure interaction represented by human vocal folds in airflow

The paper deals with the human vocal folds vibration excited by the fluid flow. The vocal fold is modelled as an elastic body assuming small displacements and therefore linear elasticity theory is used. The viscous incompressible fluid flow is considered. For purpose of numerical solution the arbitrary LagrangianEuler method (ALE) is used. The whole problem is solved by the finite element method (FEM) based solver. Results of numerical experiments with different boundary conditions are presented.


Introduction
The fluid-structure interaction is a complex and very interesting problem with many, not only technical applications.Similar physical description of interaction between fluid flows and elastic bodies can be used in various studies dealing for example with airfoil or bridge stability as well as in investigation of blood flow in compliant arteries or in biomechanics of human voice production.This paper is concentrated on the biomechanics of human vocal folds, see e.g.[1].Human voice is one of basic human being's characteristics and it plays an important role in the quality of a human life.The air flow from lungs excites the vibrations of the human vocal folds which in reverse influence air flow and vice versa.This problem was studied in many papers, see e.g.[2], [3].Therefore we speak about coupled problem.The human voice originates in basic sound produced by the vocal fold vibration.The attention of this paper is focused on fluid-structure interaction (FSI) problem and comparison of different boundary conditions.First, the interaction problem is described in the terms of linear elasticity theory and the incompressible viscous fluid flow modelled by the Navier-Stokes equations.Then the mathematical model is derived with the help of the arbitrary Lagrangian-Eulerian (ALE) method which enables to handle with the time-dependent domain.The numerical algorithm is based on the finite element method (FEM).For the fluid solver the P1-bubble/P1 elements and for the elastic part piecewise linear P1 elements are implemented.For the coupled problem the fully implicit scheme is used.This implemented algorithm is verified on few test cases when different boundary conditions are studied.

Mathematical model
For the sake of simplicity our study is restricted to a 2D model problem.Scheme of the model domain is shown in a e-mail: jan.valasek1@fs.cvut.czFlow model.Just for this moment we consider firmly given fluid domain Ω f independent on time.Then the motion of the viscous incompressible Newtonian fluid in Eulerian coordinates is described by the Navier-Stokes equations where v denotes the fluid velocity and p is the kinematic pressure, ν f is the kinematic viscosity of the fluid and g f is the force vector.ALE method.During the interaction the vocal fold changes it's shape, thus the fluid domain is also deformed and the ALE method is suitable approach to resolve it.It is the generalization of the Eulerian and the Lagrangian description.
The basic assumption of this method is the existence of diffeomorphism A t : These conditions ensure that mapping A t maps deformed fluid domain in time t back to reference shape.It follows thereof the similarity with Lagrangian method, but on contrary ALE method does not describe the real motion of the particles.
In addition mapping A t must meet the conditions This method allows us to define velocity of deformation of the fluid domain w D as where ŵD is quantity defined on Ω f re f as ŵD (X, t) = ∂ ∂t A t (X), t ∈ (0, T), X ∈ Ω f re f .Afterwards we can establish so called ALE derivative where w D (x, t) is given by (3).For more information see e.g.[4].
Then the substitution of the ALE derivative into Eg.(1) leads to Navier-Stokes equations in ALE form, see e.g.[5] (5) The formulation of the fluid problem is completed by the initial conditions and boundary conditions, specially for solid walls and the interface for the inlet we consider two possibilities: where x ∈ Γ f In , t ∈ (0, T) and n f denotes outer unit normal to ∂Ω f t .And finally boundary condition at the outlet is where x ∈ Γ f Out , t ∈ (0, T) and p re f is a reference pressure.Only one version of each condition by Eqs. ( 9) and ( 10) will be chosen for simulations.
Elastic structure.Vocal folds are modelled as an 2D compliant structure with the assumption of small displacements.It's deformation is given by, see e.g.[6] where u = (u 1 , u 2 ) is the displacement vector, ρ s is the structure density, the tensor τ s i j is the Cauchy stress tensor, the vector f s is the volume density of an acting force and X = (X 1 , X 2 ) are the reference coordinates.Using the generalized Hook law for an isotropic body, the stress tensor τ s i j can be expressed as where I is the unit matrix and e s jk is the strain tensor for small displacements e s jk = 1 2 and λ s , μ s are Lame coefficients given by the Young modulus of elasticity E s and the Poisson's ratio σ s , see [6].
Eq. ( 11) is supplied with the following initial and boundary conditions is the unit outer normal to ∂Ω s re f and q s (X, t) = (q s 1 , q s 2 ) is the aerodynamic force vector described later.
Coupling conditions.The fluid-structure interaction is a coupled problem.It means that solutions in the fluid and the structure domain are connected via boundary conditions.Moreover, the location of the interface Γ W t at time t is not known and it depends on establishing force equilibrium between aerodynamic and elastic forces, which are dependent on the position and the velocity of the interface.
The shape of the interface Γ W t in time t is the result of force equilibrium between fluid and structure.It's position is given by the deformation u, also So called dynamic condition expresses the aerodynamic force q s acting on the structure interface where ∂x i ) for i, j ∈ {1, 2} are components of the fluid stress tensor and n f (x) = (n f 1 , n f 2 ) denotes the components of the unit normal (here to the interface Γ W t ) oriented out of Ω f t .On the other hand the kinematic condition given by Eq. ( 8) enforces the correspondence between fluid motion and elastic body deformation on boundary Γ W t .

Numerical approximation
The presented mathematical model represented by the partial differential equations ( 11), ( 5) is discretized in space by the finite element method and in time by the Newmark and the backward differentiation formula of second order (BDF2), respectively.The time interval [0, T] is divided into the equidistant partition with a constant time steps Δt.
Elastic structure.The solution u is sought in the space Dir }, and W k,p (Ω) denotes the Sobolev's space.The equation ( 11) can be reformulated in the weak sense where (•, The solution u of Eq. ( 17) must satisfy that equation for all Φ ∈ V.For further details see [7].Furthermore after approximation of the space V by finite dimensional space V h ⊂ V the discrete solution u h (x, t) can be written as a linear combination of basis functions Φ j (x) ∈ V h which yields where α(t) = (α j )(t) denotes unknown time dependent coefficient of linear combination of basis function Φ j and the elements of the matrices K = (k i j ), M = (m i j ) and vector b(t) = (b i )(t) are given by . The Lagrange finite elements of the first order are used based on a triangulation, which give the first order of accuracy in space.The Newmark's method is used for numerical solution of ordinary differential equations (19), which results for a proper parameter choice to the second order accurancy in time, for application details see [7].

ALE method.
In every time step the ALE mapping must be constructed.For that purpose the pseudoelastic solver is used.It is elastic body deformation solver with unphysical Lame parameters, which ensures better distribution of mesh deformation in the vicinity of very deformed part of the mesh, e.g.around the interface.See [8].
The discrete version of mesh deformation is dealt with BDF2.
Flow model.Firstly we discretize the equations (5) in time and then in space.For time discretization we used BDF2 scheme suitable for ordinary differential equations of first order.So the ALE derivative is approximated by where for the fixed time instant t n+1 we denote In next we omit time index and we shall write Ω f := Ω f t n+1 .The application of the approximation in Eq. ( 5) leads to the scheme Then we continue with the weak formulation, where the velocity solution v in time t n+1 is sought in the functional space W 1,2 (Ω f ) and q ∈ M = L 2 (Ω f ).Further, the space of the test function X = X × X depends on choice of boundary condition.For Eq. (9 a)) the space is defined by If alternative condition Eq. (9 b)) is chosen then test functions can be nonzero on Γ f In .In the next it is written for better arrangement v n+1 as v.The weak formulation in space is acquired by multiplication of the first equation ( 21) by ϕ ∈ X and integration over the whole domain Ω f , that results in the following equations With help of Green's theorem third and fourth term can be modified The first terms of right side of Eq. (23) will be used in weak formulation and the sum of second ones form the boundary condition where the definition of test space and the knowledge that pressure is exact up to arbitrary constant was used.This condition can be prescribed on Γ f Out or also on Γ f In if the simulation driven by pressure gradient is requested.
The weak formulation then reads For derivation of condition (10 b)) we divide term where the boundary term was split to positive and negative part: (•) + = max(•, 0), (•) − = − min(•, 0), together with w D = 0 on ∂Ω f \Γ W t and ϕ = 0 on Γ W t .This explains an additional term in Eq. (10 b) ) and changes Eq. ( 25) This formulation suppresses backward intake of flow on the outlet part of the boundary and the term 1 2 ((v increases stability of final scheme.For details see [9].The FEM then approximates spaces X and M by the finite dimension spaces X h and M h , so the solution v ≈ v h can be expressed again as linear combination of basis functions leading from Eq. (25) to the system where β, γ are vectors of linear combination coefficients and The system of equations ( 28) is non-linear.For its solution the linearization v * h = v n is used.and the mathematical library UMFPACK is employed, see [10].In this article P1-bubble/P1 elements are used, which according to [11] satisfy the well-known Babuška-Brezzi condition.
Coupled problem.The aerodynamic forces are evaluated from the velocity and pressure values at adjacent triangles with the help of numerical quadrature.The algorithm solving the coupled problem is implemented in fully implicit form, sometimes called strongly coupled.It means that new time step is solved only after previous time step converges, see [7].

Numerical results
The numerical results for fluid flow interacting with the vocal fold model M5 suggested by paper [12] are presented.The model M5 together with the triangulation is shown in figure 2. Here, only one half of the channel was used as the computational domain with symmetric boundary conditions prescribed at the top of the fluid domain (y = 0.003 m).
Flow in the fixed channel without interaction.The flow solver was tested for flow prescribed by pressure gra-  is given as and it can be calculated with use of numerical quadrature.It's time development and Fourier analysis is shown in fffigure 4 and 5.
It is interesting that even without interaction the motion of fluid demonstrates periodic behaviour, here with dominant frequency 42 Hz.It is the frequency of the boundary layer separation and the formation of vortices.This phenomena is more obvious if the simulation is driven by pressure gradient because then the inlet velocity is not constant but also oscillatory.We also controlled the total energy E s of vibrating structure which consists of kinetic and potential energy -E s = E s kin + E s pot .These energies are approximated in discrete form as Figure 7 shows the time development of E s and figure 8 it's Fourier analysis.There are two peaks at the frequencies 54.54 Hz and 278.80 Hz.The first one corresponds to the first eigenfrequency of the vocal folds, see [7].The second one is excited just in case of higher inlet flow velocities and it agrees with seventh eigenfrequency (285.6 Hz).

Conclusion
The paper presented the formulation of the FSI problem focused on the example of human fold vibration in airflow.The formulation is based on FEM and the ALE method and the developed numerical scheme was implemented in an own program.Special attention was devoted to the type of the boundary conditions used.When the fluid flow was driven by the pressure gradient the simulated interaction showed the vocal fold vibration with nearly the same frequencies as for the inlet Dirichlet boundary condition.In addition the considered non-linear version of the do-nothing condition gives more robust solution especially suitable for problems with higher flow velocities at the outlet of the channel.
The financial support for the present project was partly provided by the Czech Science Foundation under the Grant No. 101/11/0207 and project SGS 13/174/OHK2/3T/12.

figure 1 .
figure 1. Quantities connected with fluid and structure are labelled by indeces f and s , respectively.The boundary of fluid domain ∂Ω f t at arbitrary time t consists of the inlet part Γ f In , the outlet part Γ f Out and Γ f Dir ∪ Γ W t .The part Γ f Dir represents fixed walls of the glottal channel and the time dependent boundary Γ W t representing the interface on the vibrating vocal folds.The undeformed (reference) shape of the domain is marked by index re f .The fluid flow is described within the Eulerian coordinates and later solved using the arbitrary Lagrangian Eulerian coordinates, see [4].The Lagrange coordinates are used for solution of the structure motion.

Fig. 1 .
Fig. 1.Schema of vocal folds model with boundaries of the glottal channel.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences,

Fig. 3 .
Fig. 3.The magnitude of velocity pattern around the fixed model at four different time instants t = 0.05 + j • 0.075 s, ( j = 0, 1, 2, 3) and detail at 0.35 s are shown.Figures are ordered from the left to the right.At the same time 25 pressure contour lines in range from −2.6 to 2.0 Pa m 3 /kg are depicted.

Fig. 4 .
Fig. 4. The kinetic energy of fluid in 10 −4 [J] in dependence on time in [s].

Fig. 5 .
Fig. 5. Fourier transformation of periodic section of E f kin , bottom axis is in [Hz].

Fig. 6 .
Fig. 6.Flow velocity magnitude in time 0.447 s a few steps before solver failure.Maximum of speed 15 m/s is out of color scale and 20 pressure isolines are spread between −20 and 30 Pa m 3 /kg.

Fig. 7 .
Fig. 7. Time development of E s , time is in [s] and energy in [10 −4 J].

Fig. 8 .
Fig. 8. Fourier analysis of E s with two dominant frequencies 54.55 Hz and 278.80 Hz.

Fig. 10 .
Fig. 10.Time development of total energy E s .

Fig. 11 .
Fig. 11.Fourier analysis of total energy E s with dominant frequency 52.0 Hz.