Grant No. GA19-04477S of Czech Science Foundation and by Grant No. SGS19/154/OHK2/3T/12 of the CTU in Pragu

Horáček, J.

Sváček, P.

Valášek, J.

English

University of West Bohemia Digital Library

The influence of penalization inlet boundary condition on the stability boundary

DSpace at University of West Bohemia

COMPUTATIONALCHANICSME November 4 - 6, 2019Srníconference with international participation35th2019The influence of penalization inlet boundary conditionon the stability boundaryJ. Vala´sˇeka, P. Sva´cˇeka, J. Hora´cˇekbaFaculty of Mechanical Engineering, Czech Technical University in Prague, Karlovo na´m. 13, Praha 2, 121 35bInstitute of Thermomechanics, Czech Academy of Sciences, Dolejsˇkova 5, Praha 8, 182 00The flow-induced vibrations of human vocal folds (VFs) are interesting and complex phe-nomenon with number of possible practical applications, see e.g. [2, 3]. Mathematically itrepresents the fluid-structure interaction (FSI) problem which in this specific case poses manydifficulties for numerical realization. One of the major difficulty is need to consider the flowdomain time evolution especially during the closing phase, see [5], followed by the questionof a suitable inlet boundary condition, see [4] etc. We present here the two-dimensional FSImodel, see Fig. 1, where linear elastic problem is coupled to the incompressible Navier-Stokesequations in the arbitrary Lagrangian-Eulerian (ALE) form. The special attention is paid to theflow inlet boundary conditions (BCs), where the penalization boundary condition is introducedas a promising alternative to the usually used Dirichlet boundary condition or to the prescribedpressure drop between inlet and outlet.Fig. 1. The scheme of FSI domain composed of fluid domain Ωft and the elastic body domain ΩsStructure. The motion of structure with density ρs is given by partial differential equationsρs∂2ui∂t2− ∂τsij∂xj= 0 in Ωs × (0, T ), (1)where u(x, t) = (u1, u2) is sought displacement and τ sij are the components of the Cauchy stresstensor [1]. The small displacements are assumed and the elastic body is modelled as isotropic.Fluid flow. In order to incorporate the effects of the fluid computational domain, the ALEmethod is used. The motion of the viscous incompressible fluid in a time dependent domain Ωftis modelled by the Navier-Stokes equations written in the ALE formDAvDt+ ((v −wD) · ∇)v − νf∆v +∇p = 0, div v = 0 in Ωft , (2)212where v(x, t) denotes the fluid velocity, p is the kinematic pressure and νf is the kinematic fluidviscosity, further the term DADtdenotes so called ALE derivative and the term wD is domaindeformation velocity, see e.g. [1].The system of equations (2) is equipped with the zero initial condition and appropriate BC,see [6]. Especially, on the inlet ΓfIn the three different inlet BC are considereda) v(x, t) = vDir(x, t) for x ∈ ΓfIn,b) (p(x, t)− pin)~nf − νf ∂v∂~nf(x, t) = −12v(v · ~nf )− for x ∈ ΓfIn, (3)c) p(x, t)~nf − νf ∂v∂~nf(x, t) = −12v(v · ~nf )− + 1(v − vDir) for x ∈ ΓfIn,where the vector nf = (nfj ) denotes the outward unit normal to the boundary ∂Ωf and α+ =max{0, α}, α− = min{0, α}, see [6]. The condition (3 c) is weakly imposed Dirichlet boundarycondition, which fulfilment is enforced with the aid of penalization coefficient .Numerical model. The FSI problem given by Eqs. (1) and (2) together with appropriateinitial and boundary conditions are discretized in space by the finite element method and in timeby the Newmark method. The partitioned approach with strong coupling is chosen. Especially,the modified streamline-upwind/Petrov-Galerkin method is used for stabilization of flow solver.For implementation details, see [6].Numerical results. Numerical results of flow-induced vibration of vocal folds (VFs) withthe full channel configuration are presented. The VF geometry is based on articles [4] and [6].The constant time step ∆t is chosen as 2.5 × 10−5 s, the densities are set to ρs = 1000 kg/m3and ρf = 1.185 kg/m3, the kinematic fluid viscosity is νf = 1.47 × 10−5 m2/s. The elasticparameters are chosen as follows: Young modulus Es = 8 kPa and Poisson ratio σ = 0.4.Comparison of different inlet boundary conditions. The following four simulations withdifferent inlet boundary conditions are compared, namely:a) Case VEL with prescribed Dirichlet BC (3 a) with vdir = (2.1, 0.0) m/s,b) Case PEN-S with penalization BC (3 c), where velocity vdir is imposed by  = 10−5,c) Case PEN-W with penalization BC (3 c), where velocity vdir is imposed by  = 5×10−4,d) Case PRES with given pressure drop ∆p = pin − pref = 400 Pa in (3 b).The inlet flow velocity and the pressure drop are displayed in Fig. 2. While the inlet velocity isfor Dirichlet BC constant and for the prescribed pressure drop highly oscillating, in cases PEN-S and PEN-W we see the moderately oscillation behaviour of inlet velocity. The amplitude ofoscillation increases with the higher values of penalization parameter  (1is decreasing).The pressure drop in the case PRES is constant, see Fig. 2 (right), whereas in cases VEL,PEN-S and PEN-W it oscillates with exponentially increasing amplitude. This can be expectedbecause the inlet velocity vDir exceeds critical velocity of flutter instability vflutter. Then the VFoscillation amplitude gradually grows and the airflow pressure increases as the channel cross-section becomes smaller. The increase of pressure drop is most rapid for case VEL, the casesPEN-S and PEN-W are time delayed. The simulations in all four cases end by the fluid flowsolver failure caused by too distorted fluid computational mesh near the top of the VFs.The gap denoting the distance between vocal folds is plotted in Fig. 3. The graduallyclosing gap corresponds well with pressure drop behaviour in cases VEL, PEN-S and PEN-W.In the case PRES the VF oscillation develops despite constant pressure drop nevertheless thedevelopment of large VF oscillation took longer time. We found explanation in the fact, thatprescribed ∆p lays more closer to stability boundary contrary to given inlet velocity vDir, whichsubstantially exceeds the stability boundary.213Fig. 2. The numerically simulated inlet airflow velocity and pressure difference between inlet and outletof the channel are shown on the left and on the right, respectively. The simulation cases VEL, PEN-S,PEN-W and PRES are compared.Fig. 3. Time development of the gap numerically simulated for cases VEL, PEN-S, PEN-W and PRES(The graph envelope is undulated due to too low sampling rate for saving the data for drawings.)Comparison of simulations with hemi-larynx and full larynx configuration. The simula-tions with the full larynx geometry are compared with the hemi-larynx geometry as performedin [6], where the flow symmetry along x-axis is assumed in order to considerably reduce com-putational time. The numerical approximation of FSI problem for both configurations at onetime instant are shown in Fig. 4.Fig. 4. The airflow velocity magnitude shown in similar time instant during closing phase of VF oscil-lation cycle. The simulation case PEN-W (full larynx configuration) is shown left, the simulation withhemi-larynx configuration is also computed with prescribed penalization BC (3 c) enforced with the aidof  = 5× 104.214Finally, the critical velocity of flutter instability vflutter is determined by the successivelyincreasing the prescribed airflow inlet velocity until the unstable VF vibration regime occurs.Further, we studied the dependence of the determined critical velocity on the penalizationparameter, see Fig. 5. The dependency of velocity vflutter on the parameter  in the range1 × 10−10 <  < 5 × 104 is compared for both configurations, i.e., the geometry of full lar-ynx and hemi-larynx. The both dependencies show the similar behaviour however exact valuesslightly differs. Nevertheless the computationally cheaper simulation of the hemi-larynx con-figuration can be used to estimate quite well the critical velocity of the full larynx configuration.Fig. 5. The dependencies of critical velocity of flutter instability on the penalization parameter. The de-pendency of simulations with the full larynx configuration is plotted by blue curve, while the dependencyof simulations with the hemi-larynx configuration is denoted as vcritical,hemi.AcknowledgmentsAuthors are grateful for support provided by Grant No. GA19-04477S of Czech Science Foun-dation and by Grant No. SGS19/154/OHK2/3T/12 of the CTU in Prague.References[1] Feistauer, M., Sva´cˇek, P., Hora´cˇek, J., Numerical simulation of fluid-structure interaction prob-lems with applications to flow in vocal folds, In: T. Bodna´r, G. P. Galdi, S. Necˇasova´ (eds.),Fluid-structure interaction and biomedical applications, Birkhauser, 2014, pp. 312-393.[2] Hora´cˇek, J., Radolf, V., Laukkanen, A. M., Experimental and computational modeling of theeffects of voice therapy using tubes, Journal of Speech, Language, and Hearing Research 62 (7)(2019) 2227-2244.[3] Mittal, R., Zheng, X., Bhardwaj, R., Seo, J. H., Xue, Q., Bielamowicz, S., Toward a simulation-based tool for the treatment of vocal fold paralysis, Frontiers in Physiology 2 (19) (2011) 1-15.[4] Sva´cˇek, P., Hora´cˇek, J., Finite element approximation of flow induced vibrations of human vocalfolds model: Effects of inflow boundary conditions and the length of subglottal and supraglottalchannel on phonation onset, Applied Mathematics and Computation 319 (2018) 178-194.[5] Sva´cˇek, P., Vala´sˇek, J., Mathematical modelling and numerical simulation of flow induced vibra-tions of vocal folds model with collisions, AIP Conference Proceedings 2116 (2019) No. 030003,doi: 10.1063/1.5113987.[6] Vala´sˇek, J., Sva´cˇek, P., Hora´cˇek, J., On suitable inlet boundary conditions for fluid-structure inter-action problems in a channel, Applications of Mathematics 64 (2) (2019) 225-251.215

http://dspace5.zcu.cz/bitstream/11025/36360/1/Val%c3%a1%c5%a1ek.pdf

The influence of penalization inlet boundary condition on the stability boundary

Abstract

Similar works

Full text

Available Versions

University of West Bohemia Digital Library

DSpace at University of West Bohemia