Acta Polytechnica https://doi.org/10.14311/AP.2021.61.0155 Acta Polytechnica 61(SI):155–162, 2021 © 2021 The Author(s). Licensed under a CC-BY 4.0 licence Published by the Czech Technical University in Prague NUMERICAL SOLUTION OF FLUID-STRUCTURE INTERACTION PROBLEMS WITH CONSIDERING OF CONTACTS Petr Sváček Czech Technical University in Prague, Faculty of Mechanical Engineering, Department of Technical Mathematics, Karlovo nám. 13, 121 35 Praha 2, Czech Republic correspondence: petr.svacek@fs.cvut.cz Abstract. This paper is interested in the mathematical modelling of the voice production process. The main attention is on the possible closure of the glottis, which is included in the model with the concept of a fictitious porous media and using the Hertz impact force The time dependent computational domain is treated with the aid of the Arbitrary Lagrangian-Eulerian method and the fluid motion is described by the incompressible Navier-Stokes equations coupled to structural dynamics. In order to overcome the instability caused by the dominating convection due to high Reynolds numbers, stabilization procedures are applied and numerically analyzed for a simplified problem. The possible distortion of the computational mesh is considered. Numerical results are shown. Keywords: Aeroelasticity, finite element method, incompressible fluid. 1. Introduction The voice production mechanism is a complex pro- cess consisting of a fluid-structure-acoustic interac- tion problem, where the coupling between fluid flow, viscoelastic tissue deformation and acoustics is cru- cial, see [1]. The so-called phonation onset (flutter instability) for a certain airflow rate with a certain prephonatory position leads to the vocal folds oscilla- tion. The important aspect of the phenomena is the glottis closure (glottis is the narrowest part between the vibrating vocal folds). The considered problem can be mathematically de- scribed as a problem of fluid-structure interaction with the involvement of the (periodical) contact problem of the vocal folds. In order to include the interactions of the fluid flow with solid body deformation an the contact problem, a simplified model problem is consid- ered. This model is similar to the simplified two mass model of the vocal folds of [2], see also the aeroelastic model in [3]. In this paper the mathematical model is introduced and the numerical approximation of the problem is described, where the residual based stabilization is used for the incompressible flow model. The simplified lumped vocal fold model is considered, where the contact is treated with the aid of the Hertz impact forces. In the flow model, the contact is considered with the combination of a suitable modification of the inlet boundary condition and the concept of a fictitious porous media flow. Attention is paid to the details of the finite element approximations with the aid of the ALE conservative method. The solution of the system of equations is discussed, with attention on the projection method and on the discrete projection method. Numerical tests are presented. 2. Mathematical model In this section a simplified model fluid-structure prob- lem is considered. The domain occupied by fluid Ωt is shown in Figure 1. The fluid flow is coupled with the elastic structure deformation of a simplified two degrees of freedom model of a vocal fold. 2.1. Flow problem First, the air flow is modelled by the system of the Navier-Stokes equations (cf. [4]) written in the ALE conservative form (cf. [5]) 1 JA DA(JAρu) Dt + ρ∇· ((u−wD) ⊗u) = divτf, ∇·u = 0, (1) where u = (v1,v2) is the fluid velocity vector, ρ is the constant fluid density, τf = (τfij) is the fluid stress ten- sor given by τf = −pI + 2µD, the symmetric gradient tensor D = (dij) is given by D(u) = 12 (∇u + ∇ Tu), p denotes the pressure and µ > 0 is the constant fluid viscosity. Further, by At the Arbitrary Lagrangian- Eulerian mapping is denoted which maps at any time t ∈ [0,T] the reference configuration Ωref = Ω0 onto the current configuration Ωt, by JA the Jacobian of this mapping is denoted, wD denotes the domain velocity, and D Au Dt is the ALE derivative, i.e. the derivative with respect to the reference configuration Ωref. For the system (1) the initial and boundary con- ditions are prescribed. The boundary conditions are prescribed on the boundary ∂Ωft of the compu- tational domain formed by mutually disjoint parts ∂Ωft = ΓI ∪ ΓS ∪ ΓO ∪ ΓWt, where ΓI denotes the inlet, ΓO the outlet, ΓS the axis of symmetry and ΓWt denotes either the fixed or the deformable wall. The standard boundary conditions are prescribed at 155 https://doi.org/10.14311/AP.2021.61.0155 https://creativecommons.org/licenses/by/4.0/ https://www.cvut.cz/en Petr Sváček Acta Polytechnica 0 L L L 2 I O g (t ) Γ ΓΓ Ω Wt t Figure 1. The computational domain Ωt with speci- fication of the boundary parts. 3 2 1 1 w 2 w 1 2 2 1 k k m F m m F CoG Figure 2. Two degrees of freedom model (with masses m1, m2, m3) in displaced position (displace- ments w1 and w2) The acting aerodynamic forces F1 and F2 are shown. Γt = ΓWt∪ΓWf formed of the fixed wall ΓWf (where wD = 0) or the moving wall ΓWt - i.e. the surface of the vocal fold model and at the axis of symmetry ΓS a) u = wD on ΓWt, b) u2 = 0,−τ f 12 = 0 on ΓS, (2) as ΓS is chosen to be located at line y = const (in computations we choose y = 0). Furthermore, at the inlet ΓI and at the outlet part ΓO the boundary conditions are formally prescribed as c) ρ 2 (u ·n)−u−n ·τf = 1 ε (u−uI) on ΓI, d) ρ 2 (u ·n)−u−n ·τf = prefn on ΓO, (3) where n denotes the unit outward normal vector to ∂Ωft , uI is a prescribed inlet velocity, pref is a refer- ence pressure value (pref = 0 in what follows), ε > 0 is a penalization parameter and α− denotes the negative part of a real number α. Here, the boundary condition (3c) weakly imposes the Dirichlet boundary condition u = uI with the aid of a penalization parameter ε. 2.2. Structure model The motion of the vocal fold is modelled as a motion of a rigid body governed by the displacements w1(t) and w2(t) of the two masses m1 and m2, respectively (see Fig. 2). The equation of motion (see [3] for details) reads Mẅ + Bẇ + Kw = −F , (4) where M is the mass matrix of the system, K is the diagonal stiffness matrix of the system characterized by spring constants k1,k2, and B = ε1M + ε2K, is the matrix of the proportional structural damping, ε1, ε2 are the constants of the proportional damping. The mass matrix is given by M = ( m1 + m34 m3 4 m3 4 m2 + m3 4 ) , K = ( k1 0 0 k2 ) , (5) where m1,m2,m3 are the masses shown in Fig. 2. The components of F = (F1,F2)T are the aerodynamical forces (downward positive) acting at masses m1 and m2, see Fig. 2. 2.3. Coupling conditions The aerodynamical forces F1,F2 are computed with the aid of the aerodynamical lift force L(t) and the aerodynamical torsional moment M(t) acting on the surface of the structure ΓWt. The aerodynamical lift force and the aerodynamical torsional moment are evaluated with the aid of the mean (kinematic) pressure p and the mean flow velocity u = (u1,u2) as the integrals over the surface of the airfoil L = −l ∫ ΓW t τ f 2jnj dS, M = l ∫ ΓW t τ f ijnjr ort i dS, (6) where l denotes the depth of the profile section, and the vector rort has components rort1 = −(x2 −xEA2 ), and rort2 = x1 −xEA1 , with (xEA1 ,xEA2 ) being the posi- tion of the structure elastic axis. The displacement of the structure surface ΓWt is determined in terms of w1,w2, and it is used as the boundary condition for the displacement of any point of the computational domain Ωref0 onto the domain Ωt. Moreover, the domain velocity at ΓWt is determined using the w1(t) and w2(t) obtained from solution of the ordinary differential equations (4). 2.4. Treatment of the contact problem In order to close of the fluid computational domain, the idea of a fictitious porous media flow is employed. It means that during the closing phase of the vocal folds, the part of the structure domain (denoted by Ωpt ) is assumed to be a domain of porous media flow, see Fig. 3. The flow in the domain Ωpt is then assumed to be governed by the modified equations ρ JA DA(JAu) Dt + ρ∇·((u−wD)⊗u) + αu = divτf, (7) where the coefficient α corresponds to the artificial porosity of the medium, see [6]. The coefficient in the porous media flow is usually chosen as α = µ P where µ is the dynamical viscosity of the fluid and P is an artificial permeability coefficient, see [6]. In practical computations equation (7) is solved in the whole computational domain Ωt with α being zero at all points outside of the porous media domain ΩPt and being a positive constant (in this paper α was chosen as 107) in the interior of ΩPt . 156 vol. 61 Special Issue/2021 Numerical solution of FSI problems with contacts Ω Γ t Wt Figure 3. The detail of the porous media flow domain Ωpt . The porous media domain is determined by the following procedure illustrated in Figure 3. First, the actual gap g(t) at time t is computed as distance of ΓWt from the axis of symmetry ΓS and compared to a prescribed minimal gap threshold gmin. If g(t) ≥ gmin, then the glottis is open and no porous media is used (i.e. ΩPt is empty and α ≡ 0). For g(t) < gmin the y-displacement of the points of the surface of ΓWt is modified in order to not violate the condition g(t) ≥ gmin, i.e. these points are vertically displaced by such a shift −∆w, which makes the actual gap equal to g(t) = gmin. For the structural model the Hertz impact forces are used as specified in [3]. In this case, the right hand side F of Eq. (4) is modified by the addition of the Hertz impact force FH. The Hertz impact force is then distributed to the both components of F depending on the position xmax of the impact. The model of Herz impact forces is given as FH = kHδ3/2(1 + bHδ̇), kH ≈ 4 3 E 1 −µ2H √ r, where δ stands for the penetration of the vocal fold through the contact plane, bH is a damping factor (here set to zero), r is the radius of the osculating circle (i.e. inverse of the curvature), E is Young’s modulus of elasticity of the vocal fold and µH is its Poisson’s ratio, for values see [3]. 2.5. Mesh deformation First, due to the motion of the structure, the trans- formation of the computational domain (and mesh as well) needs to be constructes at every time instant tk. Here it is represented by the construction of the ALE mapping At at every time instant t = tk. The ALE mapping At is sought in the form of displacements d = d(ξ,t) of points ξ of the original configuration, i.e. At(ξ) = ξ + d(ξ,t). The boundary condition for d is then known for every point of the boundary ∂Ω: at ∂Ω \ ΓWt the displacement d = 0 and at ΓWt it is determined by the already known (or estimated) position of the structure surface ΓWt in terms of the displacements w1 and w2. In order to determine d in Ω a boundary value problem is solved, see e.g. [7]. 3. Numerical approximation In this section let us consider the computational do- main Ωt and the ALE mapping At to be already known and At to be sufficiently smooth at every time instant t ∈ I = (0,T). A similar assumption is made about the domain velocity wD. 3.1. Weak formulation In order to introduce the function spaces for velocity and pressure, we start with defining the spaces for velocity and pressure at the reference configuration. We shall seek the velocity-pressure pair U = (u,p) at any time instant t in the space Wt = Vt ×Qt, where Qt = L2(Ωt) and Vt = {ϕ ∈ H1(Ωt) : ϕ ·n = 0 at ΓS}. Further, the weak formulation is derived using the standard approach, i.e. we multiply equations (1) by test functions z and q, integrate over Ωt and use Green theorem. However, in the finite element con- text such test functions are usually time independent, whereas here the test functions are defined over the time dependent domain Ωt and naturally must be chosen as time dependent. To this end we first define the test functions on the reference domain Ωref0 , i.e. we denote X ref = {ϕ ∈ H1(Ω ref 0 ) : ϕ = 0 at Γ W ref 0 , ϕ ·n = 0 at ΓS }. (8) Next, these (reference) test functions are defined on the current configuration Ωt with the aid of the ALE mapping At, i.e. the test function z = z(x,t) can be defined using a reference test function zref ∈ X ref by z(x,t) = zref (Atn+1 (A −1 t (x))) ∀x ∈ Ωt, t ∈ (0,T). (9) We shall denote the space of such test functions by X , i.e. X = {z : z satisfy (9) for some zref ∈ X ref}. (10) The space X is for any t subspace of Vt, and the test functions from X are time independent when transformed on the reference domain. This means that their ALE derivative equals zero, i.e. for z ∈X we have DAtz Dt = 0. In what follows by the symbol (·, ·)M the dot prod- uct in L2(M) or L2(M) is denoted. The weak form is then derived in the standard form: we take test 157 Petr Sváček Acta Polytechnica function V = (z,q) ∈ X ×Qt multiply the first equa- tion (1) by test functions z and the second by q, sum together, integrate over Ωt and use Green’s theorem for viscous terms and the pressure gradient. Thus we arrive to the weak form: Find U = (u,p) ∈ Wt such that for any t ∈ (0,T) u satisfy the boundary condition (2a) and d dt (ρu,z)Ωt + c(U; U,V ) + B(U; U,V ) (11) b(U,V ) + b∗(U,V ) + d(U,V ) = LΓ(V ), holds for any V = (z,q) ∈ X ×Qt. The forms in Eq. (11) are defined for any U = (u,p) ∈ Wt, U = (u,p) ∈ Wt and V = (z,q) ∈ X ×Qt as follows: the form d in Eq. (11) is defined by b(U,V ) = (∇·u,q)Ωt, b ∗(U,V ) = −(∇·z,p)Ωt, d(U,V ) = ∫ Ωt 2µdij(u)dij(z) dx the boundary forms b,LΓ are given as B(U; U,V ) = 1 ε (u,z)ΓI + (ρ 2 (u ·n)+u,z ) ΓI∪ΓO LΓ(V ) = 1 ε (uI,z)ΓI − ∫ ΓO pref (n ·z)dS and the convective term is given by the skew- symmetric trilinear form c (here we abbreviate w = u−wD) c(U; U,V ) = ( ρ 2 ((w ·∇)u,z ) Ωt − ( ρ 2 ((w ·∇)z,u ) Ωt − ( ρ 2 (∇·wD )u,z ) Ωt . The time derivative terms arises from the identity∫ Ωt 1 JA DA(JAu) Dt ·zdx = d dt (u,z)Ωt , (12) which follows from (see [8]) DAJA Dt = JA(divwD). (13) In what follows an equivalent integral form of equation (13) shall be used in the form d dt (∫ Vt dx ) = ∫ ∂Vt wD ·ndS (14) for any volume Vt whose motion is fully determined by the ALE mapping, i.e. Vt = At(V ref 0 ). 3.2. Time discretization For the purpose of time discretization an equidistant partition tj = j∆t of the time interval I with the constant time step ∆t > 0 is considered and we denote the approximations of velocity and pressure by uj ≈ u(·, tj) and pj ≈ p(·, tj) for j = 0, 1, . . . , uj ∈ Vtj and pj ∈. Similarly, by JAj and w j D the Jacobian of the ALE mapping At at time instant tj and the domain velocity at time instant tj is denoted. In what follows we shall describe the discretization at a fixed time instant tn+1. The time derivative on the right hand side of equation (12) is approximated at t = tn+1 formally by the second order backward difference formula 1 JA DA(JAu) Dt |tn+1 ≈ (15) 1 JAn+1 3JAn+1un+1 − 4JAn un + JAn−1un−1 2∆t or more precisely the time derivative term from (11) is approximated by d dt (∫ Ωt ρu ·zdx ) |tn+1 ≈ M(U,V ) −LM (V ) (16) for U = (u,p) := (un+1,pn+1), V = (z,q) and with M(U,V ) = 3 2∆t ( un+1,z ) Ωtn+1 , (17) LM (V ) = 2 ∆t (un,z)Ωtn − 1 2∆t ( un−1,z ) Ωtn−1 (18) where the test function z ∈ X is a time dependent function defined by (9) with a given steady reference function. Thus, the time discretized weak formulation reads: Find U = (u,p) := (un+1,pn+1) such that u satisfy the boundary condition (2a) and a(U; U,V ) = L(V ) (19) holds for any V = (z,q) ∈X ×Qtn+1, where a(U; U,V ) = M(U,V ) + c(U; U,V ) + B(U; U,V ) + b(U,V ) + b∗(U,V ) + d(U,V ), (20) and L(V ) = LM (V ) + LΓ(V ). 3.3. Finite element approximations The spaces Vtn+1 and X are approximated using the finite element subspaces Vh and X h constructed over an admissible triangulation T∆ of the domain Ω, re- spectively. Similarly, the pressure space Qtn+1 is ap- proximated by its finite element subspace Qh con- structed again over T∆. Here, the Taylor-Hood finite elements are used, i.e. the spaces of continuous piece- wise quadratic velocities and continuous piecewise linear pressures are used, i.e. the velocities are sought in Vh = {ϕ ∈ C(Ω) : ϕ ∈ P 2(K)∀K ∈T∆}∩ V. the space of the test functions is given by X h = {ϕ ∈ C(Ω) : ϕ ∈ P 2(K)∀K ∈T∆}∩ X . and the trial/test pressure space is given as Qh = {ϕ ∈ C(Ω) : ϕ ∈ P 1(K)∀K ∈T∆}. 158 vol. 61 Special Issue/2021 Numerical solution of FSI problems with contacts The finite element approximations of u and p are then sought in the finite element spaces Vh = Xh×Qh constructed over an admissible triangulation τh of the computational domain Ωft : Find an approximate so- lution Uh = (u,p) ∈ Wh such that Eq. (19) holds for any test function Vh = (z,q) ∈ X h×Qh. Furthermore, this formulation is stabilized using the SUPG/PSPG stabilization terms together with the div-div stabiliza- tion terms given as S(U; U,V ) = ∑ K∈T∆ δK ( ρ 3u 2∆t −µ4u + ρ (w ·∇) u + ∇p, Φ(U; V ) ) K F(U; V ) = ∑ K∈T∆ δK ( ρ 4ũn − ũn−1 2∆t , Φ(U; V ) ) K P(U,V ) = ∑ K∈T∆ τK ( ∇·u,∇·z ) K , where the modified test function is given by Φ(U; V ) = ((u−wD) ·∇)z +∇q and δK, τK are suitably chosen stabilization parameters, see e.g. The stabilized discrete formulation then reads: Find U = (u,p) ∈ Wh such that aS(U; U,V ) = LS(U; V ) (21) holds for any test function V = (z,q) ∈ X h ×Qh, where aS(U; U,V ) = a(U; U,V ) + P(U,V ) + S(U; U,V ), and LS(U; V ) = L(V ) + F(U; V ). 3.4. Linearization In order to solve the nonlinear problem (21), the sequence of the linearized problem is solved until it converges to a sufficient precision. We start with an approximation U0 and for k = 0, 1, . . . solve the linearized problems: Find U =: Uk+1 such that aS(Uk; U,V ) = LS(Uk; V ) (22) holds for any V = (z,q) ∈ X h ×Qh. Using the finite element base functions the discretization leads to the system of linear equations in the form of a saddle point problem, i.e. in the form( A B∗ BT Q )( αu αp ) = ( f g ) , (23) where the matrix A = M + D + C + AB + AS consists of the mass matrix M, the diffusion matrix D, the convection matrix C, the boundary terms matrix AB and the stabilization matrix AS part arising from the forms M, d, c, B and S, respectively. Matrix B∗ is the discrete gradient operator. It corresponds to the form b∗ including the stabilization terms, BT corresponds to the discrete divergence operator realized by the form b plus stabilization. The matrix Q follows fully from the stabilization S, i.e. in particular for the case with no stabilization the matrix Q = 0. The solution of such a problem is difficult, see e.g. [9]. This is caused by the presence of the continuity equation, which can be treated at the continuous level by the approach based on the projection method. On the other hand, this can be understood as an approximation of the solution of the system (23), or more precisely as preconditioned iterative solution of this system. 4. Solution of the discrete problem 4.1. Projection method The projection method is based on the Helmholtz- Hodge decomposition of any vector field v, i.e. v = vdiv + ∇Ψ, where vdiv is the divergence free field, see [10]. In this section, for the sake of simplicity, we shall restrict ourselves to the case of the fixed computational domain Ω and to the case of the first order discretization in time. This means that here we shall discuss the solution of the Navier-Stokes equations in the form ρ ∂u ∂t + ρ(u ·∇)u + ∇p−µ4u = 0, ∇·u = 0, where the time derivative term is replaced by the backward difference formula ρ ∂u ∂t ≈ ρ un+1 −un ∆t . Here, the problem (24) can be sought using the segre- gated approach, see [11]. This can be formally written as I. Solve momentum equations for ũ ũ−un ∆t + (ũ ·∇)ũ−ν4ũ = f, (24) II.Project ũ on the space of (discrete) divergence free functions: Solve pressure equation ∇· (∇p) = ∇· ( ũ ∆t ) , (25) III.Update the velocity field by adding the pressure gradient: un+1 − ũ ∆t + ∇p = 0. (26) Let us mention that the steps I. and II. need to be equipped with suitable boundary conditions. Due to splitting the coupled equations, the boundary con- ditions as presented in (3-2) require a modification. Particularly the boundary conditions (2) contain pres- sure, which is not available in the step I, but can be used from the previous time step. The pressure equation needs to be equipped by the Neumann type boundary condition, where the compatibility condi- tion needs to be satisfied for the existence of a solution, 159 Petr Sváček Acta Polytechnica which is unique to a constant. Nevertheless, due to splitting these two steps, the splitting error arises, for overview see e.g. [12]. In the considered problem, the main difficulty is in the realization of the non-standard boundary conditions (2,3). To this end we shall con- sider another possibility based on the preconditioned solution of the discrete problem (23). 4.2. Discrete projection methods and preconditioning Let us consider the system of linear equations written in the form Mu = b. where the matrix M is given as M = ( A B BT 0 ) , b = ( f g ) . (27) This is the well studied case, where the matrix A is (possibly) a non-symmetric positive definite and the matrix B has full rank due to the Babuška-Brezzi inf-sup condition satisfied. The matrix M can be factorized as M = ( I 0 BT A−1 I )( A 0 0 S )( I A−1B 0 I ) (28) where S is the pressure Schur complement given as S = −BT A−1B. (29) The inverse matrix to the matrix M given (28) is able to compute the inverse matrix to the matrix A and the inverse to the Schur complement matrix S. In order to solve the problem (27), the three sub-steps can be followed similarly as in Section 4.1. I. Solve the momentum equations for the intermediate velocity field and update the right hand side for pressure by subtracting the “divergence” of the intermediate velocity field, i.e. ũ := A−1f, g̃ := g −BTũ. II.Solve pressure equation Sp = g̃ III.Update the velocity field by adding a pressure gradient component to the momentum equations u = ũ−A−1(Bp) In this process, the solution of the system with ma- trix A and the matrix S needs to be realized. Here, matrix A is a mass matrix perturbed with convection- diffusion. The stabilization terms and the solution of the system with the matrix A can be realized effec- tively. However the solution of the system with the Schur complement matrix S must be realized itera- tively. This is why the above approach is used only as a preconditioner for GMRES method, where the matrix S is replaced by a suitable approximation, see e.g. [13]. Figure 4. Detail of the mesh for the almost closed glottal part. In the considered discretization, the matrix of the system is given by M = ( A BT −B D ) , b = ( f g ) . (30) where the matrix A is (possibly) a non-symmetric positive definite, the matrix B has full rank and matrix D is positive semidefinite. In this case the matrix M can be factorized as M = ( I 0 −BA−1 I )( A 0 0 S1 )( I A−1BT 0 I ) (31) where S1 is the pressure Schur complement given as S1 = S + BA−1BT . 5. Numerical Results First, the Oseen problem is considered −ν4u + b ·∇u + ∇p + αu = f, ∇·u = 0 in the computational domain Ω = (0, 1)2. The prob- lem is equipped with the Dirichlet boundary condition u = b prescribed at ∂Ω. Here, α = 0 is used, b = (sin(πx),−πy cos(πx)) and the right hand side f is chosen in such a way that u = (sin(πx),−πy cos(πx)) is solution of the Oseen problem. The computations were performed for different values of the viscosity coefficient nu. First, the convergence of the Galerkin finite element approximations uGh to the exact solu- tion u = b is investigated, p(x,y) = sin(πx) cos(πy) for ν = 0.05 (here, relatively high viscosity was cho- sen in order to obtain stable Galerkin approximations even on coarser meshes). For an approximation of the flow problem, the Taylor Hood finite elements were used. The errors in the H1 norm are shown in Table 1, where by hmax the length of the maximal edge in the triangulation is denoted. These results are compared to the results of the stabilized formulation of the same problem, which shows that the used residual based stabilization does not pollute the solution, see Table 2. The convergence orders in both cases agree well with the theoretical estimate. For the stabilized method such convergence rates are well preserved for the val- ues ν = 10−3, . . . , 10−6 with a slow down observed only for coarse grid configuration. 160 vol. 61 Special Issue/2021 Numerical solution of FSI problems with contacts hmax H 1(u) H1(v) H1(p) 0.333174 0.148971 0.278814 0.824015 0.166358 0.0294769 0.0389521 0.306831 0.0881204 0.00751673 0.00970303 0.155735 0.0449673 0.00178417 0.00225841 0.0781441 0.0230627 0.000444434 0.00055994 0.038375 0.0118955 0.000107972 0.000139096 0.0192135 Order 2.14 2.17 1.07 Table 1. Convergence of Galerkin FE method to the solution of the Oseen problem with convergence order estimate hmax H 1(u) H1(v) H1(p) 0.333174 0.148971 0.278814 0.824015 0.158053 0.0385103 0.0479786 0.332202 0.0877183 0.00912475 0.0113478 0.162359 0.0451748 0.00239609 0.00279304 0.0821471 0.0235271 0.000593443 0.000695859 0.0402822 0.0119951 0.000148245 0.000174687 0.0201917 Order 2.06 2.05 1.03 Table 2. Convergence of Stabilized FE method to the solution of the Oseen problem with convergence order estimate hmax H 1(u) H1(v) H1(p) 0.333174 0.0895357 0.149882 0.764274 0.158053 0.0211092 0.033792 0.322219 0.0877183 0.00602738 0.00931535 0.160369 0.0451748 0.00158538 0.00240926 0.0810965 0.0235271 0.000394708 0.000603108 0.0397932 Order 1.94 2 1.16 Table 3. Convergence of Stabilized FE method to the solution of the Navier-Stokes problem with con- vergence order estimate -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 t[s] -0.05 -0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04 0.05 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 t[s] Figure 5. The aeroelastic response in terms of w1(in mm, top) and w2 (in mm, bottom) of the structure for flow velocity U∞ = 0.65 m/s - phonation onset -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 1.9 1.95 2 2.05 2.1 2.15 2.2 2.25 2.3 t[s] -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 1.9 1.95 2 2.05 2.1 2.15 2.2 2.25 2.3 t[s] Figure 6. The aeroelastic response in terms of w1(in mm, top) and w2 (in mm, bottom) of the structure for flow velocity U∞ = 0.65 m/s - phonation phase with almost periodical oscillations -0.0006 -0.0004 -0.0002 0 0.0002 0.0004 0.0006 0 0.05 0.1 0.15 0.2 t[s] -0.0004 -0.0003 -0.0002 -0.0001 0 0.0001 0.0002 0.0003 0.0004 0 0.05 0.1 0.15 0.2 t[s] Figure 7. The aeroelastic response of the structure in terms of w1(in mm, top) and w2 (in mm, bottom) for flow velocity U∞ = 0.70 m/s - phonation onset -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t[s] -0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.2 0.25 0.3 0.35 0.4 0.45 0.5 t[s] Figure 8. The aeroelastic response in terms of w1(in mm, top) and w2 (in mm, bottom) of the structure for flow velocity U∞ = 0.70 m/s - phonation onset 161 Petr Sváček Acta Polytechnica 0.1 5.6 11.1 16.6 22.1 27.6 33.1 38.6 0.1 5.6 11.1 16.6 22.1 27.6 33.1 38.6 0.1 5.6 11.1 16.6 22.1 27.6 33.1 38.6 Figure 9. The flow velocity magnitude during the opening and closing phase for the inlet flow velocity U∞ = 0.65 m/s Let us emphasize, that as the convection b equals the exact solution u, the Dirichlet problem for Navier- Stokes equations can be formulated with the same analytical solution. This problem was solved as an Navier-Stokes problem with a known analytical solu- tion. One can really notice that b = u, and thus the given u is the analytical solution of the Navier-Stokes problem with the known Dirichlet boundary condition. Similarly as above, the numerical approximation was compared for both the stabilized and non-stabilized method, with the results confirming the expected the- oretical order of convergence in H1-norm, see Table 3 for the stabilized method. Further, the method for the solution of FSI problem with contact was realized and tested on a benchmark problem from [14]. The results for the inflow velocity U∞ = 0.65 m/s are shown in terms of displacements w1 and w2 in dependence on time in Figures 5-6. One can easily see that the system becomes aeroelastically unstable and the so called phonation onset phenomena arises. With further continuation, the vibrating vocal folds start to be influenced by their mutual contact, see the vibrations resembling a limit cycle oscillations in Figure 6. Similar behaviour is observed also for the inflow velocity U∞ = 0.7 m/s with a much faster appearance of the contact. Figure 9 then shows the details of the flow in terms of flow velocity in the glottis region. 6. Conclusion In this paper the mathematical description of a sim- plified fluid-structure interaction problem arising in the biomechanics of voice production was presented. Special attention was paid to considering the contact of the vibrating vocal folds, which was treated with the aid of the Hertz impact forces, a suitable choice of inlet boundary conditions and the use of a fictitious porous media flow description. Attention was also given to the numerical discretization of the problem, to the solution of the linearized problem and to the realization of the gap closing. The numerical tests were shown to verify the theoretical error estimates of the applied method and the numerical results of a benchmark test problem were presented. Acknowledgements This work was supported by the Czech Science Foundation under the Grant No. 19 - 07744S. References [1] I. R. Titze. The Myoelastic Aerodynamic Theory of Phonation. National Center for Voice and Speech, U.S.A., 2006. [2] K. Ishizaka, J. L. Flanagan. Synthesis of voiced sounds from a two-mass model of the vocal coords. The Bell System Technical Journal 51:1233–1268, 1972. [3] J. Horáček, P. Šidlof, J. G. Švec. Numerical simulation of self-oscillations of human vocal folds with Hertz model of impact forces. Journal of Fluids and Structures 20(6):853–869, 2005. [4] M. Feistauer. Mathematical Methods in Fluid Dynamics. Longman Scientific & Technical, Harlow, 1993. [5] T. Nomura, T. J. R. Hughes. An arbitrary Lagrangian- Eulerian finite element method for interaction of fluid and a rigid body. Computer Methods in Applied Mechanics and Engineering 95:115–138, 1992. [6] E. Burman, M. A. Fernández, S. Frei. A nitsche-based formulation for fluid-structure interactions with contact 2018. ArXiv:1808.08758. [7] M. Feistauer, J. Horáček, P. Sváček. Numerical simulation of vibrations of an airfoil induced by turbulent flow. Communications in Computational Physics 17(1):146–188, 2015. [8] M. Feistauer, M. Křížek, V. Sobotíková. An analysis of finite element variational crimes for a nonlinear elliptic problem of a nonmonotone type. East-West J Numer Math 1:267–285, 1993. [9] H. Elman, D. Silvester. Fast nonsymmetric iterations and preconditioning for Navier-Stokes equations. Tech. Rep. No. 263, Manchester, England, 1995. [10] A. Chorin. Numerical solution of the Navier-Stokes equations. Math Comput 22:745–762, 1968. [11] A. J. Chorin. A numerical method for solving incompressible viscous flow problems. Journal of Computational Physics 2(1):12–26, 1967. [12] A. Quarteroni, A. Valli. Numerical Approximation of Partial Differential Equations. Springer, Berlin, 1994. [13] A. Wathen, D. Silvester. Fast iterative solution of stabilised Stokes systems, part II: Using general block preconditioners. SIAM Journal on Numerical Analysis 31:1352–1367, 1994. [14] P. Sváček, J. Horáček. Numerical simulation of glottal flow in interaction with self oscillating vocal folds: comparison of finite element approximation with a simplified model. Communications in Computational Physics 12(3):789–806, 2012. doi:10.4208/cicp.011010.280611s. 162 http://dx.doi.org/10.4208/cicp.011010.280611s Acta Polytechnica 61(SI):155–162, 2021 1 Introduction 2 Mathematical model 2.1 Flow problem 2.2 Structure model 2.3 Coupling conditions 2.4 Treatment of the contact problem 2.5 Mesh deformation 3 Numerical approximation 3.1 Weak formulation 3.2 Time discretization 3.3 Finite element approximations 3.4 Linearization 4 Solution of the discrete problem 4.1 Projection method 4.2 Discrete projection methods and preconditioning 5 Numerical Results 6 Conclusion Acknowledgements References