Vibrations of Nonlinear Elastic Structure Excited by Compressible Flow

: This study deals with the development of an accurate, efﬁcient and robust method for the numerical solution of the interaction of compressible ﬂow and nonlinear dynamic elasticity. This problem requires the reliable solution of ﬂow in time-dependent domains and the solution of deformations of elastic bodies formed by several materials with complicated geometry depending on time. In this paper, the ﬂuid–structure interaction (FSI) problem is solved numerically by the space-time discontinuous Galerkin method (STDGM). In the case of compressible ﬂow, we use the compressible Navier–Stokes equations formulated by the arbitrary Lagrangian–Eulerian (ALE) method. The elasticity problem uses the non-stationary formulation of the dynamic system using the St. Venant–Kirchhoff and neo-Hookean models. The STDGM for the nonlinear elasticity is tested on the Hron–Turek benchmark. The main novelty of the study is the numerical simulation of the nonlinear vocal fold vibrations excited by the compressible airﬂow coming from the trachea to the simpliﬁed model of the vocal tract. The computations show that the nonlinear elasticity model of the vocal folds is needed in order to obtain substantially higher accuracy of the computed vocal folds deformation than for the linear elasticity model. Moreover, the numerical simulations showed that the differences between the two considered nonlinear material models are very small.


Introduction
The simulation of an interaction of flow and elastic bodies plays an important role in a number of areas of science, engineering and technology. However, the FSI plays an important role also in bio-medicine. Namely, the flow in veins or the heart or flow-induced vibration of human vocal folds (VFs) producing voice is intensively studied. The problems of FSI have been studied by a number of different methods in several books (e.g., [1][2][3][4][5][6][7]).
The mathematical and numerical modeling of hemodynamics uses the fact that blood is an incompressible liquid and its flow can be described by the incompressible Navier-Stokes (N-S) equations. The simulation of the airflow in the vocal tract is usually simulated with the aid of the same system for the description of the airflow in spite of the air being a compressible gas, using the fact that airflow is slow. Maximum glottal jet velocity is ca 45 m/s and computations at low Mach numbers (M < 0.1) are difficult, so the numerical modeling of flow-induced vibrations of the VFs prevails with the usage of incompressible flow; see, e.g., [8][9][10]. Nevertheless, in the solution of the compressible N-S equations, we use a robust numerical method with respect to the Mach number (e.g., [11]). This is the reason that, in our research, we prefer to use the model of compressible flow.
As for the motion, vibrations and deformation of VFs induced by the airflow in the vocal tract, several approaches have been used. For example, in [12], the unsteadiness of the flow is modeled by a prescribed periodic motion of a part of the channel wall with large amplitudes, nearly closing the channel. The numerical solution is implemented using the finite volume method and the predictor-corrector MacCormack scheme with artificial viscosity. For computation of the acoustic output, the driven VFs oscillation can be used for modeling the acoustic sources in incompressible, unsteady flow, followed by the usage of an acoustic analogy method; see, e.g., [13]. Various lumped parameter dynamic models of the VFs with several degrees of freedom are used for modeling the phonation onset and the VFs self-sustained vibrations; see, e.g., [14][15][16][17]. The fluid flow in the glottal channel modeled by the incompressible N-S equations in the ALE form, discretized by the finite element (FE) method and in combination with the dynamic lumped model of the VFs with three degrees of freedom, was recently focused on for the introduction of a new inlet boundary condition using a penalization approach. This gives reliable results related to the flutter analysis of the system and physically real values of pressure and flow velocities when the channel is nearly closing during the VFs vibration; see [18].
In [19,20], the interaction of compressible flow with a linear elastic body is solved. A survey can be found in [21].
In the numerical solution of FSI problems with compressible flow and nonlinear dynamic elasticity, it is necessary to overcome several obstacles: the numerical solution of compressible N-S equations in time-dependent domains, the solution of nonlinear elasticity problems, the realization of the interaction of these problems and the solution of large nonlinear algebraic systems.
For the numerical solution of compressible viscous flow, one of the most attractive techniques appears to be the discontinuous Galerkin method (DGM). It is used in a number of works (see, e.g., [22][23][24][25]). The situation is more complicated for the compressible flow in time-dependent domains. One possibility is to combine the DGM with the arbitrary ALE method. In [26], the ALE-DGM is applied to the solution of airfoil vibrations by the compressible flow. In [11,19], the ALE-DGM is applied to the interaction of compressible flow with the vibrations of a linear dynamic elastic body solved by the combination of the standard FE method for the space discretization and the well-known Newmark time discretization.
In this study, we are interested in the comparison of the St. Venant-Kirchhoff and neo-Hookean nonlinear elasticity models. Because of the successful solution of compressible flow by the DGM, we discretize the elasticity problems also by the discontinuous Galerkin method. Our goal is to compare these nonlinear elasticity models and their application to the linear elasticity model in a simulation of VF vibrations excited by airflow. Both compressible flow and elasticity problems are solved by the STDGM. The interaction of flow and elastic deformation is via transmission conditions on the boundary between the flow domain and the elastic body.
As follows from above, the novelty of this study is the application of the STDGM to the solution of the compressible N-S equations in the conservative ALE form in a timedependent domain coupled with two nonlinear elasticity models; the St. Venant-Kirchhoff and the neo-Hookean model. The developed method is applied to the numerical simulation of airflow in a simplified model of the human vocal tract and flow-induced VF vibrations.
In Section 2, the definition and STDGM discretization of the flow problem are formulated. Section 3 contains the formulation and discretization of elasticity problems. The transmission between the flow and elasticity coupled problems and the description of the ALE mapping construction are contained in Section 4. Further, Section 5 contains algorithmization and realization of the discrete coupled problem, including necessary details for the iterative algorithm of computing the nonlinear elasticity discrete problem. Section 6 is devoted to testing the method for the solution of the dynamic elasticity method.
Finally, Section 7 presents numerical experiments showing the robustness of the developed method by simulations of air-flow-induced vibration of the VFs in the model of the vocal tract. The results suggest the importance of nonlinear elasticity in such a study.

Compressible Flow
First, we describe the formulation and discretization of compressible flow in a timedependent domain. We proceed briefly, because it is described in several of our previous works, such as [11,21,[26][27][28].

Continuous Flow Problem
We consider compressible flow in a bounded domain Ω t ⊂ R 2 for t ∈ [0, T]. The boundary of Ω t consists of four disjointed parts: where Γ I represents the inlet, Γ O is the outlet and boundaries Γ W and Γ W t denote impermeable fixed and moving walls, respectively.
The time dependence of the domain Ω t is taken into account by using a regular oneto-one ALE mapping of the reference domain Ω 0 onto the current configuration Ω t : A t : . Then the continuity equation, the N-S equations and the energy equation can be written in the ALE form We have R s (w, ∇w) = ∑ 2 k=1 K s,k (w) ∂w ∂x k , where K s,k (w) are 4 × 4 matrices depending on w, and f s (w) = A(w)w with A(w) = D f s (w)/Dw, see [11]. We use the following notation: p-pressure, ρ-fluid density, E-total energy, v = (v 1 , v 2 )-velocity vector, θ-absolute temperature, c v > 0-specific heat at constant volume, γ > 1-Poisson adiabatic constant, µ > 0-dynamic viscosity, λ = −2µ/3second viscosity coefficients, k > 0-heat conduction coefficient, τ V ij -components of the viscous part of the stress tensor.
Equation (1) is completed by the following thermodynamical relations for pressure and absolute temperature and equipped with initial and boundary conditions with prescribed data w 0 , ρ D , v D , z D is the velocity of a moving wall and n denotes the unit outer normal.

Discrete Flow Problem
We describe the discretization, which is used in our in-house solver. We assume that Ω t is a polygonal domain for every t ∈ [0, T]. We denote by T ht a partition of the closure Ω t into a finite number of closed triangles with disjoint interiors satisfying standard properties (see [29]). We suppose that T ht is an image of T h0 under the regular mapping "t → A t ". Moreover, we assume that the ALE mapping A t is continuous and affine in Ω 0 .
By F , we denote the system of all faces of all elements K ∈ T ht . Moreover, we introduce the sets of boundary faces F B = {Γ ∈ F ; Γ ⊂ ∂Ω t }, "Dirichlet" boundary faces F D = {Γ ∈ F B ; a Dirichlet condition is prescribed on Γ} and inner faces F I = F \ F B . Each Γ ∈ F is associated with a unit normal vector n Γ to Γ. For Γ ∈ F B , the normal n Γ has the same orientation as the outer normal to ∂Ω t . For K ∈ T ht , h K denotes the average of K and h Γ denotes the length of Γ ∈ F .
For each Γ ∈ F I , there exist two neighboring elements K (L) Γ . We use the convention that K with S r ht = {v; v| K ∈ P r (K) ∀ K ∈ T ht }, where r > 0 is an integer and P r (K) denotes the space of all polynomials on K of degree ≤ r. It is possible to see that S r Thanks to properties of the expressions in the N-S equations, similarly as in [11], the following forms are derived: We set Θ = 1, Θ = 0 or Θ = −1 and get the so-called symmetric (SIPG), incomplete (IIPG) or nonsymmetric (NIPG) version, respectively, of the discretization of viscous terms. In (7) and (8), C W denotes a positive sufficiently large constant and K T is the transposed matrix to K.
In the form (9), symbols P + g (w, n) and P − g (w, n) denote the "positive" and "negative" parts of the matrix P g (w, n) = ∑ 2 s=1 (A s (w) − z s I)n s defined in the following way. By [30], this matrix is diagonalizable. It means that there exists a nonsingular matrix T = T(w, n) such that where λ i = λ i (w, n) are eigenvalues of the matrix P g . Now we define the "positive" and "negative" parts of the matrix P g by where λ + = max(λ, 0), λ − = min(λ, 0). The boundary state w B is defined on the basis of the Dirichlet boundary conditions (3) and extrapolation: In order to avoid spurious oscillations in the approximate solution in the vicinity of discontinuities or steep gradients, we apply artificial viscosity forms. They are based on the discontinuity indicator By [ρ h ], we denote the jump of the function ρ h on the boundary ∂K, and |K| denotes the area of the element K. Then, we define the discrete discontinuity indicator and the artificial viscosity forms (see [31] with parameters ν 1 , Because of the time discretization, we consider a partition 0 with integers r, q ≥ 1. Here, P q (I m ) denotes the space of all polynomials in t on I m of degree ≤ q and the space S r ht is defined in (4). For ϕ ∈ S rq hτ , we introduce the following notation: In order to bind the solution on intervals I m−1 and I m , we augment the resulting identity by the penalty expression Furthermore, we define the prolongation w hτ (t) of w hτ | I m−1 on the interval I m .
In what follows, we introduce the notation for functions a, b defined in a set ω ⊂ R 2 . Now, the space-time DG approximate solution is defined as a function w hτ ∈ S rq hτ satisfying (20) and the following relation for m = 1, . . . , M: Remark 1. In the derivation of the discrete problem, the approximate solution and the test functions are considered as elements of the space S rq hτ . In practical computations, integrals appearing in the definitions of the formsâ h ,b h , d h , J h ,Ĵ h andβ h and also the time integrals over I m are evaluated with the aid of quadrature formulas using values of the approximate solution at discrete points of intervals I m . Therefore, the space S rq hτ is finite-dimensional, and the discrete problem is equivalent with a finite algebraic system for every m = 1, . . . , M.

Dynamic Elasticity Problem
Following [32], we consider an elastic body represented by a bounded polygonal domain Ω b ⊂ R 2 . By ∂Ω b , we denote the boundary of the domain Ω b and assume that where F −T = (F −1 ) T . By P, we denote the first Piola-Kirchhoff stress tensor, which depends on the elasticity model. It is a function of u via F: P = P(F(u)) (we can refer to [32].) Piola-Kirchhoff We need to find a displacement function u : where we prescribe the following quantities: M ≥ 0 is the damping coefficient. We consider the following cases. (a) Linear elasticity: In this case the stress tensor P(F) = σ(u) depends linearly on the strain tensor e(u) = (∇u + ∇u T )/2, i.e., e ij = 1 Here, λ b and µ b are the Lamé parameters that can be expressed with the aid of the Young modulus E b and the Poisson ratio ν b : (b) Nonlinear elasticity: In the case of nonlinear models, we introduce the Green strain tensor E ∈ R 2×2 defined by with components One possibility is the model of neo-Hookean material with the stress tensor Another possibility is the use of the St.Venant-Kirchhoff model, when the second Piola-Kirchhoff stress tensor and the first Piola-Kirchhoff stress tensor are defined by

Discretization of the Elasticity Problem
In the discretization problem, we consider the displacement u and the deformation velocity y and split the basic system into two systems of first-order in time We construct a triangulation T b h of Ω b with standard properties. The approximate solution at every time instant t ∈ [0, T] will be sought in the finite-dimensional space where s > 0 is an integer. By F b h , we denote the system of all faces of all elements K ∈ T b h and distinguish their sets of boundaries, "Dirichlet", "Neumann" and inner faces: are tensors, then we define the tensor product by The DG discretization in space is formulated with the use of the following forms. Linear elasticity form: where σ(u) is defined by (29). Here, the parameter Θ is chosen as 1, 0, −1 for SIPG, IIPG, NIPG, respectively, version of the elasticity form. Nonlinear IIPG elasticity form (Θ = 0): Penalty form: Here, C b W > 0 is a sufficiently large constant.
Right-hand side form (with Θ = 0 in the case of nonlinear elasticity): In the nonlinear case, it is not clear how to define the SIPG and NIPG versions of the elasticity forms so that the form a b h is linear with respect to the test function ϕ.

STDGM for the Structural Problem
In the time interval [0, T], we consider the same partition 0 = t 0 < t 1 < . . . < t M = T and use the same notation as in Section 2.2. An approximate solution of problems (35)-(38), i.e., the approximations of the functions u, y will be sought in the space of piecewise By s and q * , we denote positive integers representing the degrees of polynomial approximations in space and time. We introduce the one-sided limits and jump of a function ϕ ∈ [S b,sq * hτ ] 2 at time t m similarly as in (19). Now, the approximate STDG solution of problem (35)-(38) is defined as a couple u hτ , y hτ ∈ S b,sq * hτ such that The initial states u h (0−), y h (0−) ∈ S b,s h are defined by The discrete nonlinear problem is solved on every time interval [t k−1 , t k ] by the Newton method (cf. e.g., [33]). Some details are contained in Section 5. The resulting linear algebraic systems in the flow and elasticity problems are solved by the direct solver UMF-PACK (see [34]) or the GMRES method with block diagonal preconditioning (see, e.g., [35]).

Transmission Conditions
On the fluid-structure boundarỹ we consider interface conditions representing the continuity of the normal stress and velocity: (a) linear elasticity: (b) nonlinear elasticity: Here is the stress tensor of the fluid, σ b ij -stress tensor of the structure, i, j = 1, 2. n(X) = (n 1 (X), n 2 (X))-the unit outer normal to the body Ω b on Γ b N at point X.

Construction of the ALE Mapping
ALE mapping A t is constructed by means of a solution of an artificial static linear elasticity problem according to [36]. We define d = (d 1 , d 2 ) in Ω 0 as a solution of the static problem where τ a ij are the components of the artificial stress tensor τ a ij = δ ij λ a divd + 2µ a e a ij (d), e a ij (d) The Lamé coefficients λ a and µ a are related to the artificial Young modulus E a and the Poisson number ν a as in (30). The boundary conditions for d are prescribed by The solution of the problem (48) and (49) provides the ALE mapping of Ω 0 onto Ω t in the form for each time instant t. In our computations, piecewise linear approximations of the function d, and thus A t , are used.

Coupling Procedure
In the solution of the complete coupled FSI problem, it is necessary to apply a suitable coupling procedure. See, e.g., [37], for a general framework. Here, we apply the following so-called strong coupling algorithm, in which we proceed successively from one time interval [t k , t k+1 ] to the next interval [t k+1 , t k+2 ].

1.
Let us assume the approximate solutions of the flow problem and the deformation of the structure u hτ,k on the time level t k are known.

Newton Method
The elasticity form a b h (u,ϕ) given by (41) is linear with respect to ϕ, but nonlinear in u. This results in systems of nonlinear algebraic equations solved by the Newton method (see [33]), which was applied in, e.g., [38,39], where incompressible flow model and conforming FE discretization were employed.
Let f : R N → R N . We seek a solution α ∈ R N such that f (α) = 0. The Newton algorithm to obtain a solution is the following: let α (0) be an initial guess of the solution and let ε > 0 be a tolerance. For i ≥ 0 proceed as follows: 1.

Important Ingredients of the Newton Method Implementation
Let ψ i , i = 1, . . . , N = dimV, be a basis of V. The solution u hτ can be expressed as In order to apply the Newton method as defined in Section 5.1, it is necessary to differentiate the form a b h (u hτ (α),ϕ) (and subsequently the tensor P) with respect to the coefficients α. For clarity, we shall denote the gradient with respect to α by ∇ α and the gradient with respect to X by ∇ X . Obviously, and By (23), P(F) = P(∇ X X + ∇ X u) Since ∇ X (X) is the constant unit matrix I, we introduce the simplified notatioñ P(∇ X u) = P(I + ∇ X u).

Derivatives in the Case of the Neo-Hookean Material
LetP =P(u hτ (α)) = (P ij ) 2 i,j=1 be the first Piola-Kirchhoff tensor of the neo-Hookean material as defined in (33). Let u hτ (α) = (u 1 , u 2 ). From (23) and (33), we get where Let u hτ (α) = (u 1 , Let us first express the derivatives of the determinant of F with respect to the coefficient α k . If 1 ≤ k ≤ N and i := k, then and for N < k ≤ 2N, i := k − N: The derivatives ofP(u hτ (α)) with respect to the coefficient α k are given as follows: If 1 ≤ k ≤ N and i := k, then where c 1 is the same as in (69), where c 1 is the same as in (69), c 2 the same as in (76) and ∂ ∂α k (det F) is expressed in (71).

Derivatives in the Case of the St. Venant-Kirchhoff Material
LetP =P(∇u h (α)) = (P ij ) 2 i,j=1 be the first Piola-Kirchhoff tensor of the St. Venant-Kirchhoff material as defined in (34). Let u h (α) = (u 1 , u 2 ). Then, we have The derivatives of P(∇ x u h (α)) with respect to the coefficient α k are given as follows: for 1 ≤ k ≤ N, i = k: (92)

Test of the STDGM for the Dynamic Elasticity
To verify the applicability of the STDGM for studying vibration problems, the benchmark CSM3 with the St. Venant-Kirchhoff model introduced by Turek and Hron in [40] is used. We consider a 2D rectangular domain representing an elastic clamped-free beam; see Figure 1. The beam loaded only by a gravity force has length l = 0.35 m and height h = 0.02 m. The goal is the computation of free vibrations of point A at the end of the beam, see Figure 1. We apply STDGM to the solution of the benchmark with the following data: u D = 0 at the clamped end, g N = 0 on the free surface of the beam, u 0 = 0, y 0 = 0, f b = (0, −ρ b g), ρ b = 1000 kg/m 3 , g = 2 m/s 2 , c b M = 0, E b = 1.4 · 10 6 Pa, ν b = 0.4. Moreover, we choose C b W = 6 · 10 6 . According to [40], the numerically simulated waveforms of point A vibration are represented by the mean value, amplitude and frequency of oscillations of the beam, see Table 1.
The computing was realized by the STDGM for both types of nonlinear models considered in the present study using linear polynomials in space and linear time approximations with successively decreasing time steps τ; see Tables 1 and 2. The results agree with the benchmark for both nonlinear theories very well and the differences between the St. Venant-Kirchhoff and neo-Hookean elasticity models are small. We can note that both the benchmark and our results give no exactly physically true amplitudes of the beam vibrations, because no positive or negative damping was included in the model, and therefore, the simulated amplitudes should correspond to the initial position of the beam.
Finally, in Table 3, we compare results for different computational meshes obtained for the St. Venant-Kirchhoff material with piecewise linear approximation in space and time for the time step τ = 0.02.

FSI Numerical Experiments Using STDGM
Here, we present numerical results of an FSI problem considering a simplified VFs model excited by airflow. The geometry of the airflow domain Ω t , which models the simplified subglottal, supraglottal spaces and a semicircle subdomain with an outlet Γ O to a surrounding atmosphere, is given in Figure 2. The boundaries Γ W of the airflow domain Ω t are considered as the impermeable hard sidewalls of the vocal tract including the vertical segments of the semicircle at the outlet. The computational domain Ω b marks the elastic VFs with the surface Γ W t , which creates an interface with the airflow domain. The VFs are fixed at the boundaries denoted by Γ b D .  For the fluid solver, the STDGM with a polynomial approximation of degree 2 in space and degree 1 in time is used. In the case of the elasticity solver, the IIPG version of the DGM with the penalization constant C W = 500 for inner faces and C W = 5000 for boundary edges is employed. The stabilization parameters ν 1 and ν 2 from (16) are set to 0.1. The time step τ is set to 1.0 · 10 −6 s. For the first 1000 time steps, the fluid flow is computed with the fixed boundary. Then, the part Γ W t of the boundary is released and we solve the FSI problem.
The elastic VFs are modeled by isotropic material of the density ρ b = 1040 kg m −3 . The VF model is formed by four subdomains with different elastic properties, as shown in Figure 4 and Table 4. Every VF contains 5118 elements.   The initial displacement and the initial deformation velocity are set to be zero. On the bottom, right and left straight parts of the boundary, we prescribe a homogeneous Dirichlet boundary condition (25) and on the curved part of the boundary, the Neumann boundary condition (26). The damping coefficient c b M is set to 1.0 s −1 . For the solution of the dynamic elasticity problem, we employ the NIPG version of the DGM, where the penalization constant is set to C b W = 4 · 10 6 . For the solution of the static elasticity problem (48), we employ the NIPG version of the DGM with the penalization constant C W = 10 3 . Then, the DG solution of the ALE discrete problem (48) is interpolated to a continuous approximation.
The strong coupling algorithm described in Section 4.3 with the prescribed tolerance 10 −5 is used. The prescribed tolerance was usually reached after two to three coupling subiterations. Figure 5 shows the airflow velocity field in the subglottal and supraglottal regions at five time instants of the VFs self-oscillation. In these time instants, different jet declination behind the channel constriction, i.e., the so-called "Coanda effect", can be observed, with the maximum flow velocity ca. 80 m/s. It corresponds to the Mach number Ma = 0.23.
Similarly as in Figure 5, we see the distribution of the pressure in Figure 6. The maximum air pressure is approximately constant in the subglottic part of the computational model. The minima of the pressure corresponds to centers of the vortices created in the supraglottal region and traveling to the outlet. Figure 7 shows the fluid pressure fluctuations in the middles of the inlet and outlet. The value vibrations are caused by the non-stationary flow behavior. In Figure 8, we can see the time dependence of the average values across the inlet and outlet. The average mean difference between the inlet and outlet pressure is around 1 kPa, which is in the range of subglottic pressure for ordinary phonation.

Comparison of the FSI Results for Linear and Two Nonlinear Elasticity Models
This section will be devoted to the analysis of nonlinear elasticity models in comparison with linear ones. We compare the linear strain tensor e and the nonlinear strain tensor E ∈ R 2×2 , defined by (32). For the linear elasticity, the stress tensor depends on the strain tensor e = (e ij ) 2 i,j=1 , and in the case of nonlinear elasticity, the stress tensor depends on E = e + E * , where E * = (E * ij ) 2 i,j=1 . The influence of the nonlinear part of the strain tensor is given by the ratio If R ≈ 1, then the nonlinear part of the strain tensor is very small and the linear elasticity model is sufficient, but if R ≈ 0, then it is necessary to use a nonlinear elasticity model.
First, we are concerned with the analysis of the neo-Hookean model. In this case, Figure 10 shows the numerical simulation of the VFs self-oscillations from the beginning of the FSI computation at 12 time instants. Figure 11 shows in detail the deformation of the VFs at two time instants for a maximal and minimal glottal gap. In Figures 10 and 11, the case R ≈ 1 is depicted by white and the case R ≈ 0 by a dark red color. The nonlinear part of the strain tensor takes effect near the VFs surface, especially in the narrowest part of the glottal channel and on the superior surface of the VFs. Therefore, to correctly capture deformations of the VFs, it is necessary to use a nonlinear elasticity model. Further, we present the results obtained for the St. Venant-Kirchhoff nonlinear elasticity model. From the comparison of Figures 10 and 12 and details in Figures 11 and 13, we see that the results for deformations of the VFs during self-oscillations are very similar. The St. Venant-Kirchhoff model shows slightly higher influence of the nonlinearity than the neo-Hookean material model of the VFs.

Discussion and Conclusions
This paper deals with the application of the space-time discontinuous Galerkin method to the numerical solution of the compressible flow in time-dependent domains, described by the compressible Navier-Stokes system, and to the nonlinear dynamic elasticity problems. This is described by the St. Venant-Kirchhoff and the neo-Hookean models. The main novelty is the numerical simulation of the fluid-structure interaction, namely, the vocal folds vibrations excited by the compressible flow. The elastic vocal folds consist of several layers with different material characteristics.
First, the applicability of the STDGM to the solution of a nonlinear dynamic elasticity is tested on a benchmark problem published in [40], where the elastic deformation of a vibrating beam is considered. Then the coupled fluid-structure problem is numerically solved. An important part of the presented study is oriented to the solution of the question whether the linear or nonlinear elasticity models of the vocal folds are more suitable. It follows from our analysis that the use of nonlinear elasticity models are more adequate than the linear model. The differences between the results obtained by the two nonlinear material models are very small.
In future studies, the identification of the generated acoustic signal and a remeshing in the case of the full glottal channel closure during the vocal folds oscillation period should be analyzed.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript:

ALE
arbitrary Lagrangian-Eulerian (method) DGM discontinuous Galerkin method FE finite element FSI fluid-structure interaction N-S Navier-Stokes STDGM space-time discontinuous Galerkin method VFs vocal folds