Next Article in Journal
Analysis of the Vibrations of Operators’ Seats in Agricultural Machinery Using Dynamic Substructuring
Next Article in Special Issue
Development of Parameters towards Voice Bifurcations
Previous Article in Journal
Potential for Farmers’ Cooperatives to Convert Coffee Husks into Biochar and Promote the Bioeconomy in the North Ecuadorian Amazon
Previous Article in Special Issue
Effects of Vertical Glottal Duct Length on Intraglottal Pressures in the Convergent Glottis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Vibrations of Nonlinear Elastic Structure Excited by Compressible Flow

1
Faculty of Nuclear Sciences and Physical Engineering, Czech Technical University in Prague, Trojanova 13, 120 00 Praha 2, Czech Republic
2
Faculty of Mathematics and Physics, Charles University, Sokolovská 83, 186 75 Praha 8, Czech Republic
3
Institute of Thermomechanics, Academy of Sciences of the Czech Republic, Dolejškova 5, 182 00 Praha 8, Czech Republic
4
Departement Technik, Ostschweizer Fachhochschule, Werdenbergstrasse 4, 9471 Buchs, Switzerland
*
Author to whom correspondence should be addressed.
Appl. Sci. 2021, 11(11), 4748; https://doi.org/10.3390/app11114748
Submission received: 11 April 2021 / Revised: 14 May 2021 / Accepted: 17 May 2021 / Published: 21 May 2021
(This article belongs to the Special Issue Computational Methods and Engineering Solutions to Voice II)

Abstract

:
This study deals with the development of an accurate, efficient and robust method for the numerical solution of the interaction of compressible flow and nonlinear dynamic elasticity. This problem requires the reliable solution of flow 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 fluid–structure interaction (FSI) problem is solved numerically by the space-time discontinuous Galerkin method (STDGM). In the case of compressible flow, 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 airflow coming from the trachea to the simplified 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.

1. 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 time-dependent 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.

2. Compressible Flow

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

2.1. 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: Ω t = Γ I Γ O Γ W Γ W t , 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 one-to-one ALE mapping of the reference domain Ω 0 onto the current configuration Ω t : A t : Ω ¯ 0 Ω ¯ t . Next, we define the domain velocity z ˜ ( X , t ) = t A t ( X ) , t 0 , T , X Ω 0 , z ( x , t ) = z ˜ ( A t 1 ( x ) , t ) , t 0 , T , x Ω t and the ALE derivative of the vector function w = w ( x , t ) , where x Ω t and t 0 , T : D A D t w ( x , t ) = w ˜ t ( X , t ) , where w ˜ ( X , t ) = w ( A t ( X ) , t ) , X Ω 0 , x = A t ( X ) . Then the continuity equation, the N-S equations and the energy equation can be written in the ALE form
D A w D t + s = 1 2 g s ( w ) x s + w div z = s = 1 2 R s ( w , w ) x s ,
where w = ( ρ , ρ v 1 , ρ v 2 , E ) T R 4 , g s ( w ) = f s ( w ) z s w , f s ( w ) = ( ρ v s , ρ v 1 v s + δ 1 s p , ρ v 2 v s + δ 2 s p , ( E + p ) v s ) T , R s ( w , w ) = ( 0 , τ s 1 V , τ s 2 V , τ s 1 V v 1 + τ s 2 V v 2 + k θ x s ) T , s = 1 , 2 ,   τ i j V = λ δ i j div v + 2 μ d i j ( v ) , d i j ( v ) = 1 2 v i x j + v j x i , i , j = 1 , 2 .
We have R s ( w , w ) = k = 1 2 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 ) / D w , 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 μ / 3 —second viscosity coefficients, k > 0 —heat conduction coefficient, τ i j V —components of the viscous part of the stress tensor.
Equation (1) is completed by the following thermodynamical relations for pressure and absolute temperature
p = ( γ 1 ) E ρ v 2 2 , θ = 1 c v E ρ v 2 2
and equipped with initial and boundary conditions
w ( x , 0 ) = w 0 ( x ) , x Ω 0 , ρ = ρ D , v = v D on Γ I , j = 1 2 i = 1 2 τ i j V n i v j + k θ n = 0 on Γ I , v = 0 on Γ W , v = z D ( t ) , θ n = 0 on Γ W t , j = 1 2 τ i j V n j = 0 , i = 1 , 2 , θ n = 0 , on Γ O ,
with prescribed data w 0 , ρ D , v D , z D is the velocity of a moving wall and n denotes the unit outer normal.

2.2. 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 h t 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 h t is an image of T h 0 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 h t . 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 h t , 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 ) , K Γ ( R ) T h t such that Γ K Γ ( R ) K Γ ( L ) . We use the convention that K Γ ( R ) lies in the direction of n Γ , and K Γ ( L ) lies in the opposite direction to n Γ . If Γ F B , then the element adjacent to Γ will be denoted by K Γ ( L ) .
Now we introduce the space of piecewise polynomial functions
S h t r = [ S h t r ] 4 ,
with S h t r = { v ; v | K P r ( K ) K T h t } , 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 h t r = { v ; v = A t ( v ^ ) , v ^ S h 0 r } . A function φ S h t r is, in general, discontinuous on interfaces Γ F I . If φ is a function defined on K Γ ( L ) K Γ ( R ) , then by φ Γ ( L ) and φ Γ ( R ) , we denote the values of φ on Γ considered from the interior of K Γ ( L ) and K Γ ( R ) , respectively, (if these values make sense) and set φ Γ = ( φ Γ ( L ) + φ Γ ( R ) ) / 2 , [ φ ] Γ = φ Γ ( L ) φ Γ ( R ) .
Thanks to properties of the expressions in the N-S equations, similarly as in [11], the following forms are derived:
a ^ h ( w ¯ h , w h , φ h , t ) = K T h t K s = 1 2 k = 1 2 K s , k ( w ¯ h ) w h x k · φ h x s d x Γ F h t I Γ s = 1 2 k = 1 2 K s , k ( w ¯ h ) w h x k ( n Γ ) s · [ φ h ] d S Γ F h t D Γ s = 1 2 k = 1 2 K s , k ( w ¯ h ) w h x k ( n Γ ) s · φ h d S Θ Γ F h t I Γ s = 1 2 k = 1 2 K k , s T ( w ¯ h ) φ h x k ( n Γ ) s · [ w h ] d S Θ Γ F h t D Γ s = 1 2 k = 1 2 K k , s T ( w ¯ h ) φ h x k ( n Γ ) s · w h d S ,
d h ( w h , φ h , t ) = K T h t K ( w h · φ h ) div z d x ,
J h ( w h , φ h , t ) = Γ F h t I Γ μ C W h Γ [ w h ] · [ φ h ] d S + Γ F h t D Γ μ C W h Γ w h · φ h d S ,
h ( w ¯ h , w B , φ h , t ) = Γ F h t D Γ μ C W h Γ w B · φ h d S Θ Γ F h t D Γ s = 1 2 k = 1 2 K k , s T ( w ¯ h ) φ h x k ( n Γ ) s · w B d S ,
b ^ h ( w ¯ h , w h , φ h , t ) = K T h t k + 1 K s = 1 2 ( A s ( w ¯ h ( x ) ) z s ( x ) ) I ) w h ( x ) ) · φ h ( x ) x s d x + Γ F h t I Γ P g + w ¯ h Γ , n Γ w h ( L ) + P g w ¯ h Γ , n Γ w h ( R ) · [ φ h ] d S + Γ F h t B Γ P g + w ¯ h Γ , n Γ w h ( L ) + P g w ¯ h Γ , n Γ w ¯ h ( R ) · φ h d S .
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 ) = s = 1 2 ( 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
P g = T Λ T 1 , Λ = diag ( λ 1 , , λ 4 ) ,
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
P g ± = T Λ ± T 1 , Λ ± = diag ( λ 1 ± , , λ 4 ± ) ,
where λ + = max ( λ , 0 ) , λ = min ( λ , 0 ) .
The boundary state w B is defined on the basis of the Dirichlet boundary conditions (3) and extrapolation:
w B = ( ρ D , ρ D v D 1 , ρ D v D 2 , c v ρ D θ Γ ( L ) + 1 2 ρ D | v D | 2 ) on Γ I ,
w B = w Γ ( L ) on Γ O ,
w B = ( ρ Γ ( L ) , ρ Γ ( L ) z D 1 , ρ Γ ( L ) z D 2 , c v ρ Γ ( L ) θ Γ ( L ) + 1 2 ρ Γ ( L ) | z D | 2 ) on Γ W t .
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
g t ( K ) = 1 h K | K | 3 / 4 K [ ρ ¯ h ] 2 d S , K T h t .
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 G t ( K ) = 0 if g t ( K ) < 1 , G t ( K ) = 1 if g t ( K ) 1 , and the artificial viscosity forms (see [31])
β ^ h ( w ¯ h , w h , φ h , t ) = ν 1 K T h t h K G t ( K ) K w h · φ h d x ,
J ^ h ( w ¯ h , w h , φ h , t ) = ν 2 Γ F h I 1 2 G t ( K Γ ( L ) ) + G t ( K Γ ( R ) ) Γ [ w h ] · [ φ h ] d S ,
with parameters ν 1 , ν 2 = O ( 1 ) .
Because of the time discretization, we consider a partition 0 = t 0 < t 1 < < t M = T of the time interval [ 0 , T ] and denote I m = ( t m 1 , t m ) , I ¯ m = [ t m 1 , t m ] , τ m = t m t m 1 , for m = 1 , , M . We define the space S h τ r q = ( S h τ r q ) 4 , where
S h τ r q = { ϕ ; ϕ ( x , t ) = i = 0 q t i ϕ i ( x ) , ϕ i S h t r , t I m , x Ω t , m = 1 , , M } ,
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 h t r is defined in (4). For φ S h τ r q , we introduce the following notation:
φ m ± = φ ( t m ± ) = lim t t m ± φ ( t ) , { φ } m = φ m + φ m .
In order to bind the solution on intervals I m 1 and I m , we augment the resulting identity by the penalty expression { w h τ } m 1 , φ h τ ( t m 1 + ) t m 1 . The initial state w h τ ( 0 ) S h 0 r is defined as the L 2 ( Ω h 0 ) -projection of w 0 on S h 0 r , i.e.,
w h τ ( 0 ) , φ h Ω t 0 = w 0 , φ h Ω t 0 φ h S h 0 r .
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
( a , b ) ω = ω a b d x ,
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 h τ r q satisfying (20) and the following relation for m = 1 , , M :
I m D A w h τ D t , φ h τ Ω t + a ^ h ( w ¯ h τ , w h τ , φ h τ , t ) d t + I m b ^ h ( w ¯ h τ , w h τ , φ h τ , t ) + J h ( w h τ , φ h τ , t ) + d h ( w h τ , φ h τ , t ) d t + I m β ^ h ( w ¯ h τ , w h τ , φ h τ , t ) + J ^ h ( w ¯ h τ , w h τ , φ h τ , t ) d t + ( { w h τ } m 1 , φ h τ ( t m 1 + ) ) Ω t m 1 = I m h ( w ¯ h τ , w B , φ h τ , t ) d t , φ h τ S h τ r q .
Remark 1.
In the derivation of the discrete problem, the approximate solution and the test functions are considered as elements of the space S h τ r q . In practical computations, integrals appearing in the definitions of the forms a ^ h , b ^ h , d h , J h , J ^ 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 h τ r q is finite-dimensional, and the discrete problem is equivalent with a finite algebraic system for every m = 1 , , M .

3. 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 Ω b = Γ D b Γ N b , where Γ D b Γ N b = . The deformation of the domain Ω b is described by the displacement u = ( u 1 , u 2 ) : Ω b × [ 0 , T ] R 2 and the deformation mapping ψ ( X , t ) = X + u ( X , t ) , X Ω b , t [ 0 , T ] . Further, we set
F = ψ , J = det F > 0 , Cof F = J ( F T ) ,
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 : Ω b × [ 0 , T ] R 2 such that
ρ b 2 u t 2 + c M b ρ b u t div P ( F ) = f b in   Ω b × [ 0 , T ] ,
u = u D in   Γ D b × [ 0 , T ] ,
P ( F ) n = g N in   Γ N b × [ 0 , T ] ,
u ( · , 0 ) = u 0 , u t ( · , 0 ) = y 0 in   Ω b ,
where we prescribe the following quantities: f b : Ω b × [ 0 , T ] R 2 is the density of the volume force, g N : Γ N b × [ 0 , T ] R 2 is the surface traction, u D : Γ D b × [ 0 , T ] R 2 is the displacement of Γ D b , u 0 : Ω b R 2 is the initial displacement, y 0 : Ω b R 2 is the initial deformation velocity, ρ b > 0 is the material density and c M b 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 i j = 1 2 u i x j + u j x i , according to the relations
t r ( e ( u ) ) = i = 1 2 u i x i ,
P ( F ) : = σ ( u ) = λ b t r ( e ( u ) ) I + 2 μ b e ( u ) .
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 = E b ν b ( 1 + ν b ) ( 1 2 ν b ) , μ b = E b 2 ( 1 + ν b ) .
(b) Nonlinear elasticity: In the case of nonlinear models, we introduce the Green strain tensor E R 2 × 2 defined by
E = 1 2 F T F I , E = ( E i j ) i , j = 1 2
with components
E i j = 1 2 u i x j + u j x i e i j linear   part + 1 2 k = 1 2 u k x i u k x j E i j nonlinear   part .
One possibility is the model of neo-Hookean material with the stress tensor
P ( F ) = μ b ( F F T ) + λ b log ( det F ) F T .
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
Σ = λ b tr ( E ) I + 2 μ b E , P ( F ) = F Σ .

3.1. 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
ρ b y t + c M b ρ b y div P ( F ) = f in   Ω b × [ 0 , T ] , u t y = 0 in   Ω b × [ 0 , T ] ,
u = u D in   Γ D b × [ 0 , T ] ,
P ( F ) n = g N in   Γ N b × [ 0 , T ] ,
u ( · , 0 ) = u 0 , y ( · , 0 ) = y 0 in   Ω b .
We construct a triangulation T h b of Ω ¯ b with standard properties. The approximate solution at every time instant t [ 0 , T ] will be sought in the finite-dimensional space
S h b , s = v L 2 ( Ω b ) ; v | K P s ( K ) , K T h b 2 ,
where s > 0 is an integer. By F h b , we denote the system of all faces of all elements K T h b and distinguish their sets of boundaries, “Dirichlet”, “Neumann” and inner faces: F h b , B = Γ F h b ; Γ Ω b , F h b , D = Γ F h b ; Γ Γ b D , F h b , N = Γ F h b ; Γ Γ N b and F h b , I = F h b \ F h b , B . For Γ F h b , the symbols n Γ , h Γ , K Γ ( L ) , K Γ ( R ) and for φ S h b , s symbols φ Γ ( L ) , φ Γ ( R ) , φ Γ and φ Γ have the same meaning as in Section 2.2.
If a = ( a i j ) i , j = 1 2 , b = ( b i j ) i , j = 1 2 are tensors, then we define the tensor product by a : b = i , j = 1 2 a i j b i j .
The DG discretization in space is formulated with the use of the following forms. Linear elasticity form:
a h b ( u , φ ) = K T h b K σ ( u ) : e ( φ ) d x Γ F h b , I Γ σ ( u ) · n · φ d S Γ F h b , D Γ σ ( u ) · n · φ d S Θ Γ F h b , I Γ σ ( φ ) · n · u d S Θ Γ F h b , D Γ σ ( φ ) · n · u d S ,
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 ):
a h b ( u , φ ) = K T h b K P F : φ d x Γ F h b , I Γ P F n · φ d S Γ F h b , D Γ P F n · φ d S .
Penalty form:
J h b ( u , φ ) = Γ F h I Γ C W b h Γ u · φ d S + Γ F h D Γ C W b h Γ u · φ d S .
Here, C W b > 0 is a sufficiently large constant.
Right-hand side form (with Θ = 0 in the case of nonlinear elasticity):
h b ( φ ) ( t ) = K T h b K f ( t ) · φ d x + Γ F h b , N Γ g N ( t ) · φ d S Θ Γ F h b , D Γ σ ( φ ) · n · u D ( t ) d S + Γ F h b , D Γ C W b h Γ u D ( t ) · φ d S .
Finally, we set A h b = a h b + J h b and
u , φ Ω b = Ω b u · φ d x .
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 h b is linear with respect to the test function φ .

3.2. 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 polynomial vector functions S h τ b , s q = [ S h τ b , s q ] 2 , where
V = S h τ b , s q = { v L 2 ( Ω b × ( 0 , T ) ) ; v | I m = i = 0 q t i φ i with φ i S h b , s , m = 1 , , M } .
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 h τ b , s q ] 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 h τ b , s q such that
I m ( ρ b y h τ t , φ h τ Ω b + c M ρ b y h τ , φ h τ Ω b + a h b ( u h τ , φ h τ ) + J h b ( u h τ , φ h τ ) ) d t + ( { y h τ } m 1 , φ h τ ( t m 1 + ) ) Ω b = I m h b ( φ h τ ) d t φ h τ S h τ b , s q ,
I m u h τ t , φ h τ Ω b y h τ , φ h τ Ω b d t + ( { u h τ } m 1 , φ h τ ( t m 1 + ) ) Ω b = 0 , φ h τ S h τ b , s q , m = 1 , , M .
The initial states u h ( 0 ) , y h ( 0 ) S h b , s are defined by
( u h ( 0 ) , φ h ) Ω b = ( u 0 , φ h ) Ω b φ h S h b , s , ( y h ( 0 ) , φ h ) Ω b = ( y 0 , φ h ) Ω b φ h S h b , s .
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 UMFPACK (see [34]) or the GMRES method with block diagonal preconditioning (see, e.g., [35]).

4. Fluid–Structure Coupling Implementation

4.1. Transmission Conditions

On the fluid–structure boundary
Γ ˜ W t = x R 2 ; x = X + u ( X , t ) , X Γ N b
we consider interface conditions representing the continuity of the normal stress and velocity:
(a) linear elasticity:
j = 1 2 σ i j b ( X ) n j ( X ) = j = 1 2 τ i j f ( x ) n j ( X ) , i = 1 , 2 , v ( x , t ) = u ( X , t ) t ,
(b) nonlinear elasticity:
P ( F ( u ( X , t ) ) ) n ( x ) = τ f ( x , t ) Cof ( F ( u ( X , t ) ) ) n ( x ) , v ( x , t ) = u ( X , t ) t .
Here, τ f = { τ i j f } i , j = 1 2 is the stress tensor of the fluid, σ i j b —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 Γ N b at point X .

4.2. 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
j = 1 2 τ i j a ( d ) X j = 0 in   Ω 0 , i = 1 , 2 ,
where τ i j a are the components of the artificial stress tensor τ i j a = δ i j λ a div d + 2 μ a e i j a ( d ) , e i j a ( d ) = 1 2 d i X j + d j X i , i = 1 , 2 . 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
d | Γ I Γ O = 0 , d | Γ W 0 \ Γ N b = 0 , d ( X , t ) = u ( X , t ) , X Γ N b .
The solution of the problem (48) and (49) provides the ALE mapping of Ω ¯ 0 onto Ω ¯ t in the form
A t ( X ) = X + d ( X , t ) , X Ω ¯ 0 ,
for each time instant t. In our computations, piecewise linear approximations of the function d , and thus A t , are used.

4.3. 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 ] .
  • 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.
  • Set u h τ , k + 1 0 : = u h τ , k , l : = 1 , and start the iterations:
    (a)
    Compute the stress tensor τ f and the aerodynamic force loading the structure and transform it to the interface Γ N b .
    (b)
    Solve the elasticity problem, compute the deformation u h τ , k + 1 l at time t k + 1 and approximate the flow domain Ω t k + 1 l .
    (c)
    Determine the ALE mapping A t k + 1 h l and approximate the domain velocity z h , k + 1 l .
    (d)
    Solve the flow problem on the approximation of Ω t k + 1 l .
    (e)
    If the variation of the displacement
    u h τ , k + 1 l u h τ , k + 1 l 1
    is larger than the prescribed difference, then set l : = l + 1 and go to (a). Otherwise, k : = k + 1 and go to (2).

5. Realization of the Discrete Nonlinear Elasticity Problem

5.1. Newton Method

The elasticity form a h b ( 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:
  • Compute the residual r ( i ) = f α ( i ) .
  • Stop iterations with α : = α ( i ) , if r ( i ) ε .
  • Compute δ α from
    α f α ( i ) δ α = r ( i ) .
  • Update α ( i + 1 ) : = α ( i ) δ α , set i : = i + 1 and go to 1.
Equation (51) is equivalent to a system of linear algebraic equations.

5.2. Important Ingredients of the Newton Method Implementation

Let ψ i , i = 1 , , N = dim V , be a basis of V . The solution u h τ can be expressed as
u h τ = u h τ ( α ) = i = 1 2 N α i ϕ i ,
where α = α i i = 1 2 N are the FE coefficients and ϕ i = ψ i , 0 for 1 i N and ϕ i = 0 , ψ i N for N < i 2 N form the basis of V 2 .
In order to apply the Newton method as defined in Section 5.1, it is necessary to differentiate the form a h b ( 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,
α k u h τ =   ψ i , 0 , 1 k N , i = k , α k u h τ =   0 , ψ i , N < k 2 N , i = k N ,
and
X u h τ = i = 1 N α i X ψ i , 0 + i = 1 N α i + N X 0 , ψ i = i = 1 N α i ψ i x 1 , i = 1 N α i ψ i x 2 i = 1 N α i + N ψ i x 1 , i = 1 N α i + N ψ i x 2 .
By (23),
P ( F ) = P ( X X + X u )
Since X ( X ) is the constant unit matrix I , we introduce the simplified notation
P ˜ ( X u ) = P ( I + X u ) .
Now, the gradient of the form a h b can be expressed as
α a h b ( u h τ ( α ) , φ ) = K T h b K α P ˜ X u h τ ( α ) : X φ d x Γ F h b , I F h b , D Γ α P ˜ X u h τ ( α ) n · φ d S + Γ F h b , I F h b , D Γ C W b h Γ α u h τ ( α ) · φ d S .
Let P ˜ X u h τ ( α ) = ( P i j ) i , j = 1 2 (here, for simplicity, we do not explicitly write the dependence of P i j on X u h τ ( α ) ) and let φ = ( φ 1 , φ 2 ) . Since
P ˜ X u h τ ( α ) : X φ = P 11 φ 1 x 1 + P 12 φ 1 x 2 + P 21 φ 2 x 1 + P 22 φ 2 x 2 ,
we have
α k P ˜ X u h τ ( α ) : X φ = α k P 11 φ 1 x 1 + α k P 12 φ 1 x 2 + α k P 21 φ 2 x 1 + α k P 22 φ 2 x 2 , α k P ˜ X u h τ ( α ) n · φ = α k P 11 n 1 + α k P 12 n 2 φ 1 + α k P 21 n 1 + α k P 22 n 2 φ 2 .
Now, for φ = ψ j , 0 , we find that
P ˜ X u h τ ( α ) : X φ = P 11 ψ j x 1 + P 12 ψ j x 2 ,
α k P ˜ X u h τ ( α ) : X φ = α k P 11 ψ j x 1 + α k P 12 ψ j x 2 ,
α k P ˜ X u h τ ( α ) n · φ = α k P 11 n 1 + α k P 12 n 2 ψ j .
Further, for φ = 0 , ψ j , we have
P ˜ X u h τ ( α ) : X φ = P 21 ψ j x 1 + P 22 ψ j x 2 ,
α k P ˜ X u h τ ( α ) : X φ = α k P 21 ψ j x 1 + α k P 22 ψ j x 2 ,
α k P ˜ X u h τ ( α ) n · φ = α k P 21 n 1 + α k P 22 n 2 ψ j .
In what follows, we express the derivatives of the tensor P ˜ .

5.3. Derivatives in the Case of the Neo-Hookean Material

Let P ˜ = P ˜ ( u h τ ( α ) ) = ( P i j ) i , j = 1 2 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
P 11 = μ b 1 + u 1 x 1 + c 1 1 + u 2 x 2 ,
P 12 = μ b u 1 x 2 c 1 u 2 x 1 ,
P 21 = μ b u 2 x 1 c 1 u 1 x 2 ,
P 22 = μ b 1 + u 2 x 2 + c 1 1 + u 1 x 1 ,
where
c 1 = λ b log ( det F ) μ b det F .
Let u h τ ( α ) = ( u 1 , u 2 ) = k = 1 2 N α k ϕ k , where ϕ k = ( ψ k , 0 ) for 1 k N and ϕ k = ( 0 , ψ k N ) for N < k 2 N .
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
α k det F = ψ i x 1 u 2 x 2 + 1 ψ i x 2 u 2 x 1 ,
and for N < k 2 N , i : = k N :
α k det F = ψ i x 2 u 1 x 1 + 1 ψ i x 1 u 1 x 2 .
The derivatives of P ˜ ( u h τ ( α ) ) with respect to the coefficient α k are given as follows: If 1 k N and i : = k , then
α k P 11 = μ b ψ i x 1 + c 2 α k det F 1 + u 2 x 2 ,
α k P 12 = μ b ψ i x 2 c 2 α k det F u 2 x 1 ,
α k P 21 = c 1 ψ i x 2 c 2 α k det F u 1 x 2 ,
α k P 22 = c 1 ψ i x 1 + c 2 α k det F 1 + u 1 x 1 ,
where c 1 is the same as in (69),
c 2 = λ b λ b log ( det F ) + μ b det F 2 ,
and α k det F is expressed in (70).
Finally, for N < k 2 N , we set i = k N and get
α k P 11 = c 1 ψ i x 2 + c 2 α k det F 1 + u 2 x 2 ,
α k P 12 = c 1 ψ i x 1 c 2 α k det F u 2 x 1 ,
α k P 21 = μ b ψ i x 1 c 2 α k det F u 1 x 2 ,
α k P 22 = μ b ψ i x 2 + c 2 α k det F 1 + u 1 x 1 ,
where c 1 is the same as in (69), c 2 the same as in (76) and α k det F is expressed in (71).

5.4. Derivatives in the Case of the St. Venant–Kirchhoff Material

Let P ˜ = P ˜ ( u h ( α ) ) = ( P i j ) i , j = 1 2 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
P 11 = μ u 1 x 2 u 2 x 1 u 2 x 2 + 1 + λ 2 u 1 x 1 + 1 u 2 x 2 + 1 2 1 + μ + λ 2 u 1 x 1 + 1 u 1 x 2 2 + u 2 x 1 2 + u 1 x 1 + 1 2 1 ,
P 12 = μ u 1 x 1 + 1 u 2 x 1 u 2 x 2 + 1 + λ 2 u 1 x 2 u 2 x 1 2 1 + μ + λ 2 u 1 x 2 u 1 x 2 2 + u 1 x 1 + 1 2 + u 2 x 2 + 1 2 1 ,
P 21 = μ u 2 x 2 + 1 u 1 x 2 u 1 x 1 + 1 + λ 2 u 2 x 1 u 1 x 2 2 1 + μ + λ 2 u 2 x 1 u 2 x 1 2 + u 1 x 1 + 1 2 + u 2 x 2 + 1 2 1 ,
P 22 = μ u 2 x 1 u 1 x 2 u 1 x 1 + 1 + λ 2 u 2 x 2 + 1 u 1 x 1 + 1 2 1 + μ + λ 2 u 2 x 2 + 1 u 1 x 2 2 + u 2 x 1 2 + u 2 x 2 + 1 2 1 .
Now let u h ( α ) = ( u 1 , u 2 ) = k = 1 2 N α k ϕ k , where ϕ k = ( ψ k , 0 ) for 1 k N and ϕ k = ( 0 , ψ k N ) for N < k 2 N . The derivatives of P ( x u h ( α ) ) with respect to the coefficient α k are given as follows: for 1 k N , i = k :
α k P 11 = μ ψ i x 2 u 2 x 1 u 2 x 2 + 1 + λ 2 ψ i x 1 u 2 x 2 + 1 2 1 + μ + λ 2 ψ i x 1 u 1 x 2 2 + u 2 x 1 2 + u 1 x 1 + 1 2 1 + 2 μ + λ 2 u 1 x 1 + 1 u 1 x 2 ψ i x 2 + u 1 x 1 + 1 ψ i x 1 ,
α k P 12 = μ ψ i x 1 u 2 x 1 u 2 x 2 + 1 + λ 2 ψ i x 2 u 2 x 1 2 1 + μ + λ 2 ψ i x 2 u 1 x 2 2 + u 1 x 1 + 1 2 + u 2 x 2 + 1 2 1 + 2 μ + λ 2 u 1 x 2 u 1 x 2 ψ i x 2 + u 1 x 1 + 1 ψ i x 1 ,
α k P 21 = μ u 2 x 2 + 1 ψ i x 2 u 1 x 1 + 1 + μ u 2 x 2 + 1 u 1 x 2 ψ i x 1 + λ u 2 x 1 u 1 x 2 ψ i x 2 + 2 μ + λ 2 u 2 x 1 u 1 x 1 + 1 ψ i x 1 ,
α k P 22 = μ u 2 x 1 ψ i x 2 u 1 x 1 + 1 + μ u 2 x 1 u 1 x 2 ψ i x 1 + λ u 2 x 2 + 1 u 1 x 1 + 1 ψ i x 1 + 2 μ + λ 2 u 2 x 2 + 1 u 1 x 2 ψ i x 2 .
For N < k 2 N , i = k N :
α k P 11 = μ u 1 x 2 ψ i x 1 u 2 x 2 + 1 + μ u 1 x 2 u 2 x 1 ψ i x 2 + λ u 1 x 1 + 1 u 2 x 2 + 1 ψ i x 2 + 2 μ + λ 2 u 1 x 1 + 1 u 2 x 1 ψ i x 1 ,
α k P 12 = μ u 1 x 1 + 1 ψ i x 1 u 2 x 2 + 1 + μ u 1 x 1 + 1 u 2 x 1 ψ i x 2 + λ u 1 x 2 u 2 x 1 ψ i x 1 + 2 μ + λ 2 u 1 x 2 u 2 x 2 + 1 ψ i x 2 ,
α k P 21 = μ ψ i x 2 u 1 x 2 u 1 x 1 + 1 + λ 2 ψ i x 1 u 1 x 2 2 1 + μ + λ 2 ψ i x 1 u 2 x 1 2 + u 1 x 1 + 1 2 + u 2 x 2 + 1 2 1 + 2 μ + λ 2 u 2 x 1 u 2 x 1 ψ i x 1 + u 2 x 2 + 1 ψ i x 2 ,
α k P 22 = μ ψ i x 1 u 1 x 2 u 1 x 1 + 1 + λ 2 ψ i x 2 u 1 x 1 + 1 2 1 + μ + λ 2 ψ i x 2 u 1 x 2 2 + u 2 x 1 2 + u 2 x 2 + 1 2 1 + 2 μ + λ 2 u 2 x 2 + 1 u 2 x 1 ψ i x 1 + u 2 x 2 + 1 ψ i x 2 .

6. 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 M b = 0 , E b = 1.4 · 10 6 Pa , ν b = 0.4 . Moreover, we choose C W b = 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 Table 1 and Table 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 .

7. 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 Γ D b .
The fluid flow problem is computed on the triangulation with 17,652 elements; see Figure 3.
Further, for the flow problem, the following input data are used:
inlet velocity v i n = 4  m s 1 ,
dynamic viscosity μ = 1.80 · 10 5  kg m 1  s 1 ,
inlet density ρ i n = 1.225  kg m 3 ,
initial outlet pressure p o u t = 97 , 611  Pa,
Reynolds number R e = ρ i n v i n H I / μ = 6941.7 ,
heat conduction coeff. κ = 2.428 · 10 2  kg m s 3  K 1 ,
specific heat c v = 721.428  m 2  s 2  K 1 ,
Poisson adiab. const. γ = 1.4 .
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 M b 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 W b = 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.
Figure 9 shows the displacement of the top of the vocal fold in horizontal ( u 1 ) and vertical ( u 2 ) directions. The numerical simulation of the VF vibrations started from zero initial conditions and static fluid forces deformed the VFs statically, predominantly in the horizontal direction. This effect shifted the mean value of VFs self-oscillations downstream. The vibration amplitudes with frequency ca. 100 Hz in the horizontal direction are dominant.

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 i j ) i , j = 1 2 , and in the case of nonlinear elasticity, the stress tensor depends on E = e + E , where E = ( E i j ) i , j = 1 2 .
The influence of the nonlinear part of the strain tensor is given by the ratio
R : = e E = e e + E .
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 Figure 10 and Figure 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 Figure 10 and Figure 12 and details in Figure 11 and Figure 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.

8. 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.

Author Contributions

Conceptualization, M.B., M.F., J.H. and A.K.; formal analysis, M.B., M.F., J.H. and A.K.; investigation, M.B., M.F., J.H. and A.K.; software, M.B. and A.K.; visualization, M.B. and A.K.; writing—original draft, M.B., M.F., J.H. and A.K.; writing—review and editing, M.B., M.F., J.H. and A.K. All authors have read and agreed to the published version of the manuscript.

Funding

The research was supported by the grant 20-01074S (M. Feistauer) and the grant 19-04477S (J. Horáček) of the Czech Science Foundation, and the research of M. Balázsová was supported by the Ministry of Education, Youth and Sports of the Czech Republic under the OP RDE grant number CZ.02.1.01/0.0/0.0/16 019/0000778/Centre for Advanced Applied Sciences, group THEORY.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the sheer large size of the dataset.

Acknowledgments

The research was supported by the grant 20-01074S (M. Feistauer) and the grant 19-04477S (J. Horáček) of the Czech Science Foundation and the research of M. Balázsová was supported by the Ministry of Education, Youth and Sports of the Czech Republic under the OP RDE grant number CZ.02.1.01/0.0/0.0/16 019/0000778/Centre for Advanced Applied Sciences, group THEORY.

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:
ALEarbitrary Lagrangian–Eulerian (method)
DGMdiscontinuous Galerkin method
FEfinite element
FSIfluid–structure interaction
N-SNavier–Stokes
STDGMspace-time discontinuous Galerkin method
VFsvocal folds

References

  1. Fung, Y.C. An Introduction to the Theory of Aeroelaticity; Dover Publications: New York, NY, USA, 1969. [Google Scholar]
  2. Dowell, E.H. Aeroelasticity of Plates and Shells; Kluwer: Dordrecht, The Netherlands, 1974. [Google Scholar]
  3. Naudasher, E.; Rockwell, D. Flow-Induced Vibrations; A. A. Balkema: Rotterdam, The Netherlands, 1994. [Google Scholar]
  4. Dowell, E.H. A modern Course in Aeroelaticity; Kluwer: Dordrecht, The Netherlands, 1995. [Google Scholar]
  5. Bisplinghoff, R.L.; Ashley, H.; Halfman, R.L. Aeroelaticity; Dover: New York, NY, USA, 1996. [Google Scholar]
  6. Paidoussis, M.P. Fluid-Structure Interactions. Slender Structures and Axial Flow, Vol I.; Academic Press: San Diego, CA, USA, 1998. [Google Scholar]
  7. Paidoussis, M.P. Fluid-Structure Interactions. Slender Structures and Axial Flow, Vol II.; Academic Press: London, UK, 2004. [Google Scholar]
  8. Mittal, R.; Erath, B.D.; Plesniak, M.W. Fluid dynamics of human phonation and speech. Annu. Rev. Fluid Mech. 2013, 45, 437–467. [Google Scholar] [CrossRef]
  9. Seo, J.H.; Mittal, R. A higher-order immersed boundary method for acoustic wave scattering and low-Mach number flow-induced sound in complex geometries. J. Comput. Phys. 2011, 230, 1000–1019. [Google Scholar] [CrossRef] [Green Version]
  10. Zheng, X.; Bielamowicz, S.; Luo, H.; Mittal, R. A computational study of the effect of false vocal folds on glottal flow and vocal fold vibration during phonation. Ann. Biomed. Eng. 2009, 37, 625–642. [Google Scholar] [CrossRef] [Green Version]
  11. Dolejší, V.; Feistauer, M. Discontinuous Galerkin Method—Analysis and Applications to Compressible Flow; Springer: Cham, Switzerland, 2015. [Google Scholar]
  12. Pořízková, P.; Kozel, K.; Horáček, J. Flows in convergent channel: Comparison of numerical results of different mathematical models. Computing 2013, 95, 573–585. [Google Scholar] [CrossRef]
  13. Šidlof, P.; Zorner, S.; Huppe, A. A hybrid approach to the computational aeroacoustics of human voice production. Biomech. Model. Mechanobiol. 2014, 14, 473–488. [Google Scholar] [CrossRef] [PubMed]
  14. Ishizaka, K.; Flanagan, J.L. Synthesis of voiced sounds from a two-mass model of the vocal cords. Bell Syst. Tech. 1972, 51, 1233–1268. [Google Scholar] [CrossRef]
  15. Liljencrants, J. A translating and rotating mass model of the vocal folds. Speech Transm. Lab. Q. Prog. Status Rep. 1991, 1, 1–18. [Google Scholar]
  16. Titze, I.R. Phonation threshold pressure: A missing link in glottal aerodynamics. J. Acoust. Soc. Am. 1992, 91, 2928–2934. [Google Scholar] [CrossRef]
  17. Story, B.H.; Titze, I.R. Voice simulation with a body cover model of the vocal folds. J. Acoust. Soc. Am. 1995, 97, 1249–1260. [Google Scholar] [CrossRef]
  18. Sváček, P.; Horáček, J. Finite element approximation of flow induced vibrations of human vocal folds model: Effects of inflow boundary conditions and the length of subglottal and supraglottal channel on phonation onset. Appl. Math. Comput. 2018, 319, 178–194. [Google Scholar] [CrossRef]
  19. Feistauer, M.; Kučera, V.; Prokopová, J. Discontinuous Galerkin solution of compressible flow in time-dependent domains. Math. Comput. Simul. 2010, 80, 1612–1623. [Google Scholar] [CrossRef]
  20. Hasnedlová, J.; Feistauer, M.; Horáček, J.; Kosík, A.; Kučera, V. Numerical simulation of fluid-structure interaction of compressible flow and elastic structure. Computing 2013, 95, 343–361. [Google Scholar] [CrossRef]
  21. Feistauer, M.; Horáček, J.; Sváček, P. Numerical Simulation of Fluid-Structure Interaction Problems with Applications to Flow in Vocal Folds. In Fluid-Structure Interaction and Biomedical Applications; Bodnár, T., Galdi, G.P., Nečasová, Š., Eds.; Springer: Basel, Switzerland, 2014. [Google Scholar]
  22. Hatrmann, R.; Houston, P. Symmetric interionr penalty DG methods for the compressible Navier-Stokes equations I: Method formulation. Int. J. Numer. Anal. Model. 2006, 1, 1–20. [Google Scholar]
  23. Dumbser, M. Arbitrary high order PNPM schemes on unstructured meshes for the compressible Navier-Stokes equations. Comput. Fluids 2010, 39, 60–76. [Google Scholar] [CrossRef]
  24. Dolejší, V. On the discontinuous Galerkin method for the numerical solution of the Navier—Stokes equations. Int. J. Numer. Methods Fluids 2004, 45, 1083–1106. [Google Scholar] [CrossRef]
  25. Dolejší, V. Semi-implicit interior penalty discontinuous Galerkin methods for viscous compressible flows. Commun. Comput. Phys. 2008, 4, 231–274. [Google Scholar]
  26. Česenek, J.; Feistauer, M.; Kosík, A. DGFEM for the analysis of airfoil vibrations induced by compressible flow. ZAMM Z. Angew. Math. Mech. 2013, 93, 387–402. [Google Scholar] [CrossRef]
  27. Feistauer, M.; Hasnedlová-Prokopová, J.; Horáček, J.; Kosík, A.; Kučera, V. DGFEM for dynamical systems describing interaction of compressible fluid and structures. J. Comput. Appl. Math. 2013, 254, 17–30. [Google Scholar] [CrossRef]
  28. Kosík, A. Fluid-Structure Interaction. Ph.D. Thesis, Department of Numerical Mathematics, Faculty of Mathematics and Physics, Charles University in Prague, Prague, Czech Republic, 2016. [Google Scholar]
  29. Ciarlet, P.G. The Finite Element Method for Elliptic Problems; North-Holland: Amsterdam, The Netherlands, 1979. [Google Scholar]
  30. Feistauer, M.; Felcman, J.; Straškraba, I. Mathematical and Computational Methods for Compressible Flow; Clarendon Press: Oxford, UK, 2003. [Google Scholar]
  31. Feistauer, M.; Kučera, V. On a robust discontinuous Galerkin technique for the solution of compressible flow. J. Comput. Phys. 2007, 224, 208–221. [Google Scholar] [CrossRef]
  32. Ciarlet, P.G. Mathematical Elasticity, Volume I, Three-Dimensional Elasticity; Volume 20 of Studies in Mathematics and Its Applications; Elsevier Science Publishers B.V.: Amsterdam, The Netherlands, 1988. [Google Scholar]
  33. Deuflhard, P. Newton Methods for Nonlinear Problems, Affine Invariance and Adaptive Algorithms; Springer Series in Computational Mathematics 35; Springer: Berlin/Heidelberg, Germany, 2004. [Google Scholar]
  34. Davis, T.A.; Duff, I.S. A combined unifrontal/multifrontal method for unsymmetric sparse matrices. ACM Trans. Math. Softw. 1999, 25, 1–19. [Google Scholar] [CrossRef]
  35. Saad, Y.; Schultz, M.H. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Statist. Comput. 1986, 7, 856–869. [Google Scholar] [CrossRef] [Green Version]
  36. Young, Z.; Mavriplis, D.J. Unstructured dynamic meshes with higher-order time integration schemes for the unsteady Navier-Stokes equations. In Proceedings of the 43rd AIAA Aerospace Sciences Meeting and Exhibit, AIAA Paper 1222, Reno, NV, USA, 10–13 January 2005. [Google Scholar]
  37. Badia, S.; Codina, R. On some fluid-structure iterative algorithms using pressure segregation methods. Application to aeroelesticity. Int. J. Numer. Methods Eng. 2007, 72, 46–71. [Google Scholar] [CrossRef]
  38. Fernandez, M.A.; Moubachir, M. A Newton method using exact jacobians for solving fluid-structure coupling. Comput. Struct. 2005, 83, 127–142. [Google Scholar] [CrossRef] [Green Version]
  39. Richter, T. Goal oriented error estimation for fluid-structure interaction problems. Comput. Methods Appl. Mech. Eng. 2012, 223–224, 28–42. [Google Scholar] [CrossRef] [Green Version]
  40. Turek, S.; Hron, J. Proposal for numerical benchmarking of fluid-structure interaction between an elastic object and laminar incompressible flow. In Fluid-Structure Interaction: Modelling, Simulation, Optimisation; Bungartz, H.J., Schäfer, M., Eds.; Springer: Berlin/Heidelberg, Germany, 2006; pp. 37–385. [Google Scholar]
Figure 1. Setup of the benchmark problem: elastic beam attached to a rigid cylinder.
Figure 1. Setup of the benchmark problem: elastic beam attached to a rigid cylinder.
Applsci 11 04748 g001
Figure 2. Computational domain at time t = 0 : L I = 20.0 mm, L g = 17.5 mm, L O = 55.0 mm, H I = 25.5 mm, H O = 2.76 mm. The radius of the semicircle subdomain is 3.0 cm.
Figure 2. Computational domain at time t = 0 : L I = 20.0 mm, L g = 17.5 mm, L O = 55.0 mm, H I = 25.5 mm, H O = 2.76 mm. The radius of the semicircle subdomain is 3.0 cm.
Applsci 11 04748 g002
Figure 3. Triangulation of the fluid domain.
Figure 3. Triangulation of the fluid domain.
Applsci 11 04748 g003
Figure 4. Nonhomogeneous model of VFs formed by four layers with triangulation. Modeled layers: 1. muscle, 2. ligament, 3. superficial lamina propria and 4. epithelium.
Figure 4. Nonhomogeneous model of VFs formed by four layers with triangulation. Modeled layers: 1. muscle, 2. ligament, 3. superficial lamina propria and 4. epithelium.
Applsci 11 04748 g004
Figure 5. Velocity field in the glottal region at five time instants ( t = 0.0074 , 0.0084 , 0.0094 , 0.0104 , 0.0114 s ) of the vocal folds self-oscillation.
Figure 5. Velocity field in the glottal region at five time instants ( t = 0.0074 , 0.0084 , 0.0094 , 0.0104 , 0.0114 s ) of the vocal folds self-oscillation.
Applsci 11 04748 g005
Figure 6. Distribution of the pressure in the glottal region at five time instants ( t = 0.0074 , 0.0084 , 0.0094 , 0.0104 , 0.0114 s ) of the vocal folds self-oscillation.
Figure 6. Distribution of the pressure in the glottal region at five time instants ( t = 0.0074 , 0.0084 , 0.0094 , 0.0104 , 0.0114 s ) of the vocal folds self-oscillation.
Applsci 11 04748 g006
Figure 7. The fluid pressure fluctuation at the middle point of inlet Γ I and at the middle point of the outlet of the vocal tract.
Figure 7. The fluid pressure fluctuation at the middle point of inlet Γ I and at the middle point of the outlet of the vocal tract.
Applsci 11 04748 g007
Figure 8. Average fluid pressure on the inlet Γ I and on the outlet Γ O .
Figure 8. Average fluid pressure on the inlet Γ I and on the outlet Γ O .
Applsci 11 04748 g008
Figure 9. Displacement of the top of the vocal fold.
Figure 9. Displacement of the top of the vocal fold.
Applsci 11 04748 g009
Figure 10. Visualization of the VFs vibrations and the ratios R of the norms of linear and nonlinear strain tensors at twelve time instants. Computed by the neo-Hookean model.
Figure 10. Visualization of the VFs vibrations and the ratios R of the norms of linear and nonlinear strain tensors at twelve time instants. Computed by the neo-Hookean model.
Applsci 11 04748 g010
Figure 11. Details of VFs deformations and the ratios R of the norms of linear and nonlinear strain tensors for the smallest and the largest glottal gap. Computed by the neo-Hookean elasticity model.
Figure 11. Details of VFs deformations and the ratios R of the norms of linear and nonlinear strain tensors for the smallest and the largest glottal gap. Computed by the neo-Hookean elasticity model.
Applsci 11 04748 g011
Figure 12. Visualization of the VFs vibrations and the ratios R of the norms of linear and nonlinear strain tensors at twelve time instants. Computed by the St. Venant–Kirchhoff elasticity model.
Figure 12. Visualization of the VFs vibrations and the ratios R of the norms of linear and nonlinear strain tensors at twelve time instants. Computed by the St. Venant–Kirchhoff elasticity model.
Applsci 11 04748 g012
Figure 13. Details of VFs deformations and the ratios R of the norms of linear and nonlinear strain tensors for the smallest and the largest glottal gap. Computed by the St. Venant–Kirchhoff elasticity model.
Figure 13. Details of VFs deformations and the ratios R of the norms of linear and nonlinear strain tensors for the smallest and the largest glottal gap. Computed by the St. Venant–Kirchhoff elasticity model.
Applsci 11 04748 g013
Table 1. Comparison of the position of point A for the reference and STDGM results with decreasing time step τ , for the St. Venant–Kirchhoff model. The displacements u 1 and u 2 are written in the format “mean value ± amplitude” and frequency. The row marked by ”ref” shows the reference benchmark values published in [40].
Table 1. Comparison of the position of point A for the reference and STDGM results with decreasing time step τ , for the St. Venant–Kirchhoff model. The displacements u 1 and u 2 are written in the format “mean value ± amplitude” and frequency. The row marked by ”ref” shows the reference benchmark values published in [40].
Method τ u 1 × 10 3 [ mm ] f [Hz] u 2 × 10 3 [ mm ] f [Hz]
ref 14.305 ± 14.305 1.0995 63.607 ± 65.160 1.0995
STDGM0.04 14.072 ± 14.043 1.0925 66.374 ± 61.499 1.0925
STDGM0.02 14.337 ± 14.316 1.0925 66.456 ± 62.556 1.0925
STDGM0.01 14.546 ± 14.526 1.0950 66.580 ± 62.994 1.0950
STDGM0.005 14.628 ± 14.608 1.0930 66.623 ± 63.153 1.0930
Table 2. Comparison of the position of point A for the STDGM with decreasing time step τ for the neo-Hookean model.
Table 2. Comparison of the position of point A for the STDGM with decreasing time step τ for the neo-Hookean model.
Method τ u 1 × 10 3 [ mm ] f [Hz] u 2 × 10 3 [ mm ] f [Hz]
STDGM0.04 14.027 ± 13.992 1.0937 66.625 ± 61.263 1.0925
STDGM0.02 14.290 ± 14.264 1.0937 66.710 ± 62.311 1.0925
STDGM0.01 14.505 ± 14.480 1.0937 66.824 ± 62.277 1.0925
STDGM0.005 14.590 ± 14.566 1.0930 66.863 ± 62.944 1.0930
Table 3. Comparison of the position of point A for STDGM with s = 1 , q = 1 for St. Venant–Kirchhoff material and different meshes defined by number of elements.
Table 3. Comparison of the position of point A for STDGM with s = 1 , q = 1 for St. Venant–Kirchhoff material and different meshes defined by number of elements.
Num. of Elem. τ u 1 × 10 3 [ mm ] f [Hz] u 2 × 10 3 [ mm ] f [Hz]
ref 14.305 ± 14.305 1.0995 63.607 ± 65.160 1.0995
7220.02 14.337 ± 14.316 1.0925 66.456 ± 62.556 1.0925
13480.02 14.117 ± 14.112 1.0962 64.508 ± 63.514 1.0962
28220.02 14.113 ± 14.110 1.0962 64.523 ± 63.518 1.0962
Table 4. Prescribed material constants for the VFs model: Young modulus, Poisson ratio and Lamé parameters for different layers, ordered from the lower layer to the upper layer. See Figure 4 for the visualization of the corresponding subdomains.
Table 4. Prescribed material constants for the VFs model: Young modulus, Poisson ratio and Lamé parameters for different layers, ordered from the lower layer to the upper layer. See Figure 4 for the visualization of the corresponding subdomains.
Layer E b ν b λ b μ b
1. layer (orange) 12 · 10 3 0.4 17,1434285
2. layer (yellow) 8 · 10 3 0.4 11,4302857
3. layer (blue) 1 · 10 3 0.495 33,110335
4. layer (red) 100 · 10 3 0.4 142,85735,714
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Balázsová, M.; Feistauer, M.; Horáček, J.; Kosík, A. Vibrations of Nonlinear Elastic Structure Excited by Compressible Flow. Appl. Sci. 2021, 11, 4748. https://doi.org/10.3390/app11114748

AMA Style

Balázsová M, Feistauer M, Horáček J, Kosík A. Vibrations of Nonlinear Elastic Structure Excited by Compressible Flow. Applied Sciences. 2021; 11(11):4748. https://doi.org/10.3390/app11114748

Chicago/Turabian Style

Balázsová, Monika, Miloslav Feistauer, Jaromír Horáček, and Adam Kosík. 2021. "Vibrations of Nonlinear Elastic Structure Excited by Compressible Flow" Applied Sciences 11, no. 11: 4748. https://doi.org/10.3390/app11114748

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop