A 2D finite-element scheme for fluid–solid–acoustic interactions and its application to human phonation

Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

A 2D finite-element scheme for fluid–solid–acoustic interactions and its application to human phonation G. Link a, M. Kaltenbacher b,*, M. Breuer c, M. Döllinger d a

Dept. of Sensor Technology, University Erlangen-Nuremberg, Paul-Gordan-Strasse 3/5, 91052 Erlangen, Germany Applied Mechatronics, University of Klagenfurt, Universitätstraße 67-69, A-9020 Klagenfurt, Austria c Dept. of Fluid Mechanics, Helmut-Schmidt University Hamburg, Holstenhofweg 85, D-22043 Hamburg, Germany d University Hospital Erlangen, Medical School, Department of Phoniatrics and Pediatric Audiology, Bohlenplatz 21, 91054 Erlangen, Germany b

a r t i c l e

i n f o

Article history: Received 20 February 2009 Received in revised form 4 June 2009 Accepted 18 June 2009 Available online 23 June 2009 Keywords: Fluid–solid interaction Aeroacoustics Fluid–solid–acoustic interaction Finite-element method Human phonation

a b s t r a c t We present a recently developed approach for the modeling of fluid–solid–acoustic interaction problems. For the efficient numerical solution of the coupled three-field problem we apply the finite-element method. The mechanical and the acoustic fields are approximated by a standard Galerkin scheme. A residual-based stabilization method is chosen for the fluid field. The interaction of the Eulerian fluid field with the Lagrangian mechanical field is based on the Arbitrary-Lagrangian–Eulerian (ALE) method and is iteratively coupled in a strong sense. The solid–acoustic interaction is based on continuum mechanics, and the fluid–acoustic coupling on Lighthill’s analogy. The new steps of our scheme are verified through validation examples. Finally, a fluid–solid–acoustic simulation of the human phonation process is presented based on a realistic model. Ó 2009 Elsevier B.V. All rights reserved.

1. Introduction In a multifield phenomenon the interactions of the physical fields involved play a crucial role. The focus of this paper lies in the continuum mechanical fields: fluid and solid mechanics as well as acoustics. Their relationship is sketched in Fig. 1. Fluid forces act thereby on a neighbouring solid, which is deformed and thus influences the velocity of the adhering fluid. Due to the deformation of the solid, the fluid domain changes and has to be adapted (grid adaption process). The fluid–acoustic interaction is described by Lighthill’s acoustic analogy and the solid–acoustic coupling by forcing continuity of surface velocities in normal direction. Fig. 1 shows the approaches chosen in the mathematical modeling of the fields and their discretization. The fluid field is modeled therein with the incompressible unsteady Navier–Stokes equations. The solid field is described by Navier’s equations, modeling linear elasticity as well as geometric nonlinearity. The acoustic sound propagation is described by the inhomogeneous wave equation based on Lighthill’s analogy and extended to account for mechanical–acoustic coupling. We restrict to such fluid–solid–acoustic interaction problems, where the influence of the acoustic field on the flow and the structure can be neglected. The aim to treat multifield problems numerically is still a growing area of research and two main challenges can be identified for * Corresponding author. E-mail address: [email protected] (M. Kaltenbacher). 0045-7825/$ - see front matter Ó 2009 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2009.06.009

the coupled field problems considered: the efficient treatment of fluid–solid interactions and the computation of the flow-induced sound (aeroacoustics). The precise and efficient computation of fluid–solid interaction problems is a highly complicated task especially for soft and/or thin-walled structures and a topic of ongoing research. Most schemes are based on the Arbitrary-Lagrangian– Eulerian (ALE) method, which allows coupling of the Eulerian fluid field with the Lagrangian mechanical field. Early works can be found in [34,37] and for current advances on this topic we refer to [68]. The main feature of ALE based approaches is to keep the topology of the fluid grid and just deform it according to the mechanical displacement along the common fluid–solid interface. To avoid any mesh deformation or even re-meshing, fixed-grid methods, especially immersed boundary (IB) methods, have recently gained great importance [44,50,54,68]. However, IB methods may suffer from the fact that there is no way to decouple the physical domain from the fictitious computational domain and/or to get rid of additional flow computation in the fictitious domain. Recently, a partitioned, iterative coupled scheme based on the eXtended FEM and DLM/FD (Distributed Lagrange Multiplier/Fictitious Domain) method has been published [69], which combines the benefits of ALE and IB methods. Towards the modeling and numerical simulation of human voice phonation, the first continuum mechanical models can be found in [2] using FEM for the computation of vocal fold vibrations. In their 2D model the mechanical field was discretized with finite-elements and the fluid forces were modeled based on the Bernoulli equation.

3322

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

Fig. 1. Fluid–solid–acoustic interaction: Arbitrary-Lagrangian–Eulerian (ALE); right-hand side (RHS); finite-element method (FEM).

A comparison between a Bernoulli and a 2D Navier–Stokes solver was discussed in [65] with a hemilarynx model (i.e. half larynx). In [63] a hemilarynx continuum mechanical model was used to clarify the causes for self-sustained vocal fold oscillations. They detected a cyclic variation of the glottis profile from a convergent to a divergent shape as the key factor for self-sustained vocal fold oscillations. In [59] a collision model with fluid–solid interaction was applied within a strongly coupled fluid–solid algorithm to treat the interactions. Full 3D coupled simulations to analyze human phonation are very challenging for current computing resources. In [31] a 3D linear elastic finite-element model of a single vocal fold with a rigid midplane surface was presented. The fluid field was eliminated and quasi-static pressure boundary conditions were applied. A 3D fluid–solid coupled model based on FEM can be found in [55]. The authors can clearly demonstrate the self-sustained vocal fold oscillations, although the grid for the flow computation was quite coarse. The first fluid–solid interaction model with fully resolved fluid can be found in [49]. They used a realistic 2D setup and applied the IB method for fluid–solid interaction. Their results are quite similar to the present ones, especially in that they were also able to obtain the so-called Coanda effect [11,17]. Therewith, the jet shows significant asymmetry and attaches stochastically to either side of the pharynx wall. This phenomenon is caused by a minor asymmetry of the jet in the geometrical configuration which induces a pressure field leading to a curvature of the jet. The curved streamlines strengthen the pressure gradient normal to the mean flow direction. Thus any initial disturbance is reinforced and the jet finally attaches to the curved wall. The second challenge for the considered three-field problem is the appropriate computation of the flow-induced sound. This topic, also known as aeroacoustics, has started with the contributions of Lighthill [46,47]. Direct solutions for acoustics and flow with the compressible Navier–Stokes equations are computationally not feasible for low Mach numbers [13,66,19]. The disparity in amplitudes, length and frequency scales, contradicts the monolithic numerical approach because the numerical grid cannot be chosen optimally for all occurring scales [66]. Hybrid methods are developed for this purpose where Lighthill’s acoustic analogy provides the theoretical background [46]. These methods are generally composed of a fluid mechanical solution step and a subsequent acoustic solution step. In general, a flow simulation is performed on a fluid domain including all relevant sound generating flow structures. Based on the results of the flow simulation, acoustic sources are computed. These sources are then used to analyze the acoustic sound propagation [5,53,20,58,42]. A comprehensive introduction to numerical methods in aeroacoustics is given in [45,66]. Towards human phonation, the aeroacoustic aspect has been first investigated in [71–73]. They described the aerodynamic generation of sound in a rigid pipe under forced vibration. The fluid–solid interaction was neglected and they focused on fluid–

acoustic coupling based on Lighthill’s acoustic analogy, which was solved with an integral method – the so-called Ffowcs Williams–Hawkins (FWH) method [22]. The results based on the FWH method were in good compliance with the results based on direct numerical simulations, which solve the compressible Navier–Stokes equations. Complementary to these studies, a theoretical approach can be found in [43]. The acoustic source model was based on a prescribed jet profile using a train of vortex rings and applied to an axisymmetric model of the vocal tract. Summarizing, we can state that many investigations considering the fluid–solid interaction on the one hand and the fluid– acoustic interaction on the other hand have been made. To the best knowledge of the authors, currently no contribution is based on the complete coupled system taking into account the fluid–solid– acoustic interaction. With the current approach both interactions are combined to achieve a numerical scheme, which is capable to simulate fluid–solid–acoustic interaction problems with special focus on human phonation. The rest of the paper is organized as follows. In Section 2, we describe the governing equations of fluid and solid mechanics as well as the field interactions considered: fluid–solid interaction, fluid– acoustic interaction and mechanical–acoustic interaction. The FE formulation of each individual partial differential equation as well as in the context of the coupled scheme is discussed in Section 3. Here we first describe in detail the used residual-based stabilized FEM for the Navier–Stokes equations. Furthermore, we discuss our mesh smoothing technique using a pseudo-structure approach with local choice of elasticity, and our aeroacoustic approach extended to account for vibrational sound. For validation we discuss in Section 4 results for a fluid–solid and a fluid–acoustic benchmark problem. In Section 5, we present the setup for the 2D simulation of the human phonation and discuss in detail the results obtained towards the flow field, mechanical vibrations of the vocal folds, and the resulting sound. A summary and conclusions are given in Section 6.

2. Physics of fluid–solid–acoustic interactions 2.1. Fluid mechanics An incompressible flow can generally be assumed for Mach numbers smaller than 0.3. Many technical and biological flows are therefore incompressible, and the governing set of partial differential equations (PDEs) is given by the momentum and the mass conservation [21]

@~ v ~ ~ v ¼ 0; þ v  rv þ rp  mD~ @t r~ v ¼ 0;

ð1Þ ð2Þ

3323

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

with p the kinematic pressure ðp ¼ P=qf Þ; m the kinematic viscosity ðm ¼ l=qf Þ and qf the density of the fluid. In order to represent moving boundaries in fluid domains, which naturally arise in fluid–solid interacting applications at the common interface Cfs , the ArbitraryLagrangian–Eulerian (ALE) approach is utilized.

interactions are coupled along the common interface, whereas fluid–acoustic interaction is a volume coupled phenomenon. For many applications it is possible to reduce the interaction to a one-way coupling. This simplification is herein assumed for fluid–acoustic and solid–acoustic interactions.

2.2. Solid mechanics

2.3.1. Fluid–solid interaction The fluid–solid interaction takes place at the domain interface Cfs and is a surface coupled interaction. At this interface kinematic and dynamic continuity have to be ensured. The complete coupled problem has to fulfill the condition that the location of the fluid– solid interface coincides for both fields

For deriving the PDE for structural mechanics, we start at Navier’s equation [38] 2

@ ~ u ~ f X þ rrs ¼ qm 2 : @t

ð3Þ

In (3) ~ f X denotes any mechanical volume force, rs the Cauchy stress u the mechanical displacetensor, qm the density of the solid and ~ ment. Now, we express rs according to Hook’s law via the linear mechanical strain lin s and the tensor of elasticity c

rs ¼ clin s

ð4Þ

~ xs;0 þ ~ u xf ¼ ~

ð5Þ

In (5) B is a differential operator, which explicitly computes for the 2D plane strain case as follows

0

@ @x

B B¼@0

@ @y

0

1

@ C @y A:

ð6Þ

@ @x

Using (4)–(6) and substituting these relations into (3) results in the PDE for linear elasticity 2

@ ~ u qm 2  BT cB~ u ¼~ f X: @t

ð7Þ

In case of large deformation, we have to extend the linear strain–displacement relation by the nonlinear terms and arrive at the Green–Lagrangian strain tensor s

s ¼

1 1 ðrX ~ u þ ðrX~ uÞT Þ þ ððrX~ uÞT rX~ uÞ; 2 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} 2

ð8Þ

lin s

where rX defines the nabla-operator with respect to the reference configuration. Furthermore, we have to transform the Cauchy stress tensor to the reference configuration resulting in the 2nd Piola Kirchhoff stress tensor T ss ¼ JF 1 d rs F d :

ð9Þ

In (9) J denotes the determinant of the deformation gradient F d , which is defined as

Fd ¼

@~ x : ~ @X

ð10Þ

Therewith, we arrive at the PDE for large deformations with respect to the reference configuration 2

@ ~ u ~ f X þ rX ðss F Td Þ ¼ qm 2 ; @t

ð11Þ

where we assume St. Venant-Kirchhoff material law to linearly relate ss and s [4]. 2.3. Field interactions The field interactions are realized through boundary conditions as well as source terms and can be generally separated into volume and surface coupled phenomena. Fluid–solid and solid–acoustic

ð12Þ

As the fluid–solid interface is impermeable and the fluid adheres to the solid, no-slip/no-penetration is assumed. The fluid–solid interaction boundary condition concerning the fluid IBVP (Initial Boundary Value Problem) is of an inhomogeneous Dirichlet type

and introduce the linear strain–displacement relation sym ~ lin u ¼ B~ u: s ¼ r

on ð0; TÞ  Cfs :

~ v¼

@~ d @t

on ð0; TÞ  Cfs ;

ð13Þ

with ~ d ¼~ u on Cfs . The movement of the fluid–solid interface leads to a change of the computational fluid domain. In order to keep the structure of the fluid grid, we will need an appropriate grid adaption, which will result in a movement of the fluid grid governed v g . By applying the ALE formulation, the conby the grid velocity ~ vective term in (1) will change to

ð~ v ~ v g Þ  r~ v ¼~ v c  r~ v:

ð14Þ

The fluid–solid boundary condition for the solid IBVP is given by an inhomogeneous Neumann condition of the form

rs  ~ n ¼ rf  ~ n;

ð15Þ

describing the equivalence of fluid stresses rf and solid stresses rs in normal direction ~ n. The fluid action on the solid according to (15) is equivalent to a force ~ f fluid . The fluid forces can be split into a pressure and a shear component

Z Z ~ pI  ~ n dC þ lðr~ v þ ðr~ v ÞT Þ  ~n dC : f fluid ¼ qf Cfs Cfs |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Pressure

ð16Þ

Shear

The fluid–solid interaction is also denoted as Dirichlet-to-Neumann problem because of the boundary condition types. 2.3.2. Fluid–acoustic interaction – Aeroacoustics Basically the compressible Navier–Stokes equations capture acoustic wave propagation phenomena in fluids. However, due to energy and length scale disparities and the preservation of the multipole character of acoustic sources, we base our approach on Lighthill’s acoustic analogy. Starting at the compressible massand momentum conservation, Lighthill’s equation in its pressure formulation is derived as [46,27]

1 @ 2 Pa  DPa ¼ r  ðr  TÞ; c2 @t 2

ð17Þ

with T the Lighthill tensor

T ¼ qf ~ v ~ v  |fflfflfflfflffl{zfflfflfflfflffl} Reynolds stress

s |{z}

Viscous stress

þ ½ðP  P0 Þ  c2 ðqf  q0 ÞI : |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

ð18Þ

Heat conduction

Thereby, P0 denotes the atmospheric pressure, c the speed of sound, qf the fluid density and q0 its atmospheric density. By considering the contribution of the three terms, simplifications to Lighthill’s tensor can be made. The Reynolds number defines the ratio

3324

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

between inertia and viscous forces. The contribution of the viscous stresses are therefore negligible for high Reynolds numbers [27,33]. In regions of ambient temperature the contribution of heat conduction is of the same order as the viscous term and can therefore be neglected as well [27]. The acoustic sources in a known flow field are finally approximated by

T  qf ~ v ~ v:

ð19Þ

In Lighthill’s acoustic analogy the acoustic decomposition is introduced only for pressure and density but not for velocity. As the mean velocity field may provide contributions to the acoustic sources, the computed acoustic pressure may have a non-zero mean component [5]. Moreover, this acoustic pressure drift may vary in space and time and has to be subtracted from the solution. 2.3.3. Solid–acoustic interaction The solid–acoustic interaction is a surface coupled phenomenon. In many technical and biological applications a vibrating solid is immersed in an acoustic fluid [41], such as e.g., electrodynamic loudspeakers, engine noises, etc. The vibrating solid causes wave propagation within the fluid and the acoustic pressure acts as a surface pressure load onto the solid. For the applications considered we can assume the solid–acoustic interaction as a one-way coupling from solid to acoustics. Along the common interface Csa the continuity of the mechanical surface v s ¼ @~u=@t and the acoustic particle velocity ~ v a have to velocity ~ meet in normal direction

~ v a  ~n ¼

@~ u ~ n on ð0; TÞ  Csa : @t

ð20Þ

In the acoustic formulation the primary unknown variable is the sound pressure P a . Therefore, the boundary condition (20) has to be formulated in terms of the pressure. The linearized Euler equation provides a relation between the sound velocity and sound pressure

@~ va ~ 1 @Pa n¼ @t q0 @~ n

ð21Þ

and we arrive at

@Pa @ 2~ u ¼ q0 2  ~ n on ð0; TÞ  Csa : @~ n @t

ð22Þ

As soon as the solid vibration is known, the acoustic boundary condition can be computed. 3. Numerical fundamentals of the fluid–solid–acoustic computational scheme The concept of semi-discretization is applied which means that spatial and temporal discretization is chosen individually and consecutively. The spatial discretization is performed by FEM, which provides a system of DAEs (Differential Algebraic Equations). The time discretization is accomplished by a FDM (Finite-DifferenceMethod) yielding a LAS (Linear Algebraic System) of equations. The time discretization of the acoustic and the solid mechanical field is performed with the Newmark method. The backward-difference (BDF2) method is applied to discretize time derivatives of the fluid field [32]. All LASs are solved with the sparse direct solver of the PARDISO library [57]. Fig. 2 displays a general setup for a fluid–solid–acoustic interaction problem and defines the computational domains, interfaces and boundaries. Therein, Xs denotes the computational domain for structural mechanics, Xf the domain for fluid mechanics and Xa ¼ Xf [ Xp for acoustics. The domain Xp is the acoustic propagation domain, where we solve a homogeneous acoustic wave equation. The fluid–solid interface is Cfs , the limiting boundary

Fig. 2. Computational setup for a fluid–solid–acoustic interaction problem.

for the fluid flow Cf and for acoustics Co . Furthermore, to simplify our derivation, we have assumed in Fig. 2 that the coupling between mechanics and acoustics just occurs along Cfs ¼ Csa (interface condition given by (22)). 3.1. Computational fluid dynamics Most flow problems are treated with the finite-volume method (FVM), which has the advantage of being mass and momentum conserving on discrete volumes. However, ensuring discrete conservation is not equivalent to small numerical errors. The velocity and pressure distributions have first priority, and hence the numerical error of these field variables is the most important property of a numerical discretization scheme. The residual-based stabilized FEM applied here eliminates the instabilities without reducing the accuracy or destroying consistency [67]. The FEM is based on a variational formulation deducible from the strong form of the incompressible Navier–Stokes equations. The formulation of the viscous term is used to obtain naturally an open domain boundary condition. Alternatively to the applied convective form v Þ, the stress-divergence form of the viscous part ð2mr  ð~ v ÞÞ ðmD~ is often used [29,30,67]. The advantage of the stress-divergence form is that the application of the integral theorem of Green to the viscous part leads to natural boundary conditions representing real physical forces. However, the physical interpretation of natural boundary conditions of the convective form implies difficulties. There are still unanswered questions concerning open boundary conditions. Especially their mathematical interpretation is not fully clarified yet and the development of open boundary conditions is an ongoing research topic [30]. The function spaces for the field variables are defined as (d denotes the space dimension) [1]

V ¼ f~ v 2 ðH1 Þd ðXf Þj~ v ¼ ~g on CfD g; Q ¼ fp 2 L2 ðXf Þg ~ of the momentum equation and for test functions w

~ 2 ðH1 Þd ðXf Þjw ~ ¼ 0 on CfD g: W ¼ fw Therewith, the weak form of (1) and (2) read as follows: Find v ; pÞ 2 ðV  Q Þ such that1 ð~

~ Þ þ ð~ ~ Þ  ðp; r  w ~ Þ þ mðr~ ~ ÞC ~ Þ  ðr  ~ ð~ v_ ; w v c  r~ v; w v ; rw v ; qÞ ¼ ~h; w fN ~ ; qÞ 2 ðW  Q Þ: 8ð w ð23Þ In order to write (23) in a more compact way, the following differential operator A is introduced 1 In the following a point over a letter denotes first-order time derivative, two points second-order time derivative.

3325

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

~ ; qgÞ ¼ ð~ ~ Þ þ ð~ ~ Þ þ mðr~ ~Þ Að~ v c ; f~ v ; pg; fw v_ ; w v c  r~ v; w v ; rw ~ Þ  ðr  ~ v ; qÞ  ðp; r  w

Table 1 for different schemes. Momentum stabilization operators Lstab M

~ ; qgÞ ¼ ð~ ~ ÞC Að~ v c ; f~ v ; pg; fw h; w fN

ð24Þ

~ ; qÞ 2 ðW  Q Þ: 8ðw

~: Lstab ¼rw C

The integral theorem of Green is applied to the viscous and the pressure term yielding a natural boundary condition containing v =@~n ¼ 0, which is known to be an appropriate open boundary @~ condition [30]. If the pressure at the boundary is assumed to be zero, which is generally possible, open boundary conditions are obtained naturally by omitting the surface integral at the right-hand side of (24). The surface integral of the pressure has to be adopted for non-zero pressure values at open boundaries. As these open boundary conditions are obtained naturally, they are also called the do-nothing boundary conditions. This formulation also enables to prescribe pressure drop conditions. A pressure drop can be established by treating the inlet and the outlet as open boundaries and by prescribing different pressure values. The open boundary can therefore be understood as an inflow–outflow condition. The standard Galerkin finite-element approach applied to the incompressible Navier–Stokes equations possesses two sources of instabilities. The first one is a result of the convective term. With increasing Reynolds numbers the impact of the convective term increases and stabilization is needed. The second instability arises from the saddle point structure of the equations and standard finite-elements do not fulfill the LBB (Ladysenskaja Babuska Brezzi) condition and therefore lead to an unstable discretization [8]. In context of FVM stabilization techniques are applied, as e.g., upwinding [21]. An enormous amount of publications exists regarding stabilized FEMs. A comprehensive overview is given in [67]. Many different stabilization methods lead to similar systems of DAEs. In order to circumvent the LBB condition and to still obtain a stable fluid mechanical solution, residual-based stabilization is applied herein for flows with moderate Reynolds numbers [9,29]. Overviews of the development can be found, e.g., in [7,16,67]. The residual-based stabilized FEM represents a group of stabilized methods including Streamline Upwind Petrov Galerkin/Pressure Stabilized Petrov Galerkin (SUPG/PSPG) [61], Galerkin Least Squares (GLS) [25] and Unusual Stabilized FEM (USFEM) [24]. In 1995, a new sight on the residual-based stabilized FEM has been initiated by introducing the variational multiscale method (VMM) [36]. The VMM provides the physical background of the additional terms appearing in the stabilized FEM [29,39]. The stabilization is based on the residual of the weak form of the incompressible Navier–Stokes equations. The idea of residual-based stabilization is to add terms including the residual of the momentum equation LM and the residual of the continuity equation LC given as

LM ¼ ~ v_ þ ~ v c  r~ v  mD~ v þ rp; v; LC ¼ r  ~

ð25aÞ ð25bÞ

which are tested with stabilization operators Lstab M;C . The stabilized variational formulation is therefore given by

~ h ; qh gÞ þ Að~ v h;c ; f~ v h ; ph g; fw

X K

~ h ÞC : ¼ ð~ hh ; w fN

~ ~ þ rq v c  rw ~  mDw ~ þ rq ~ v c  rw ~ þ mDw ~ þ rq ~ v c  rw

SUPG/PSPG [9] GLS [67] USFEM [67]

and (23) may be written as follows

sm ðLM ; Lstab M ÞK þ

X

sc ðLC ; Lstab C ÞK

K

ð26Þ

Depending on the stabilization operators Lstab M;C different methods are obtained. In Table 1, the operators of SUPG/PSPG, GLS and USFEM are listed for the momentum equation. The stabilization operator of the continuity equation denoted by grad-div stabilization was introduced in [67] and reads as follows

ð27Þ

In our scheme SUPG/PSPG together with grad-div stabilization is applied to gain a stable formulation. The crucial point of all stabilized FEM is the appropriate choice of the stabilization parameters. In (26) two stabilization parameters are included (sm and sc ) in order to be able to stabilize momentum and continuity equations individually. A review on the different choices of the stabilization parameters can be found in [67]. The appropriate choice of a parameter formulation depends on the flow problem at hand. In this paper stabilization parameters are applied, which evaluate for second-order quadrilateral elements (Q2) as

sm ¼

hK f and 2k~ v k2

sc ¼

k~ v k2 hK f; 2

ð28Þ

with k~ v k2 the L2 -norm of ~ v ; hK the discretization parameter (longest edge) and

 f¼

ReK ; 0 6 ReK < 1; 1;

ReK P 1;

:

ð29Þ

ReK denotes the element based Reynolds number defined by

ReK ¼

k~ v k 2 hK : 24m

ð30Þ

The index K denotes the elements of the grid. The inner products of the stabilization terms only have to be solved in the element interior, which is denoted by ð; ÞK in (26). In the following section the index K is omitted for simplicity. Numerical simulations of unsteady flow problems demand a time discretization. Two different philosophies have been developed to solve (26). One philosophy splits the problem into two separated sub-problems for pressure and velocity. Therewith, the overall computing effort for one iteration is smaller. However, the physical interpretation and the appropriate choice of boundary conditions in the sub-steps have not yet been clarified. Now, in the context of fluid–solid interactions, the accuracy near the walls is of particular importance. Therefore, we apply the BDF2-scheme to the unsteady flow problem (26) and obtain

    2 2 4 1 M þ DtN v nþ1 þ DtGpnþ1 ¼ M v n  v n1 : 3 3 3 3

ð31Þ

In (31) v, p denote the nodal vector of unknown velocities and pressures, n the time step, Dt the time step size and M, N, G the finiteelement matrices. The BDF2 combines a second-order accuracy with low numerical dissipation and a low numerical complexity. BDF2 is A-stable and even L-stable [32]. The time discretization leads to an algebraic system, in which the nonlinearities are preserved. Before the algebraic system can be solved, linearization needs to be applied to gain a linear algebraic system (LAS). The nonlinearity stems from the convective uÞ and its linearization in the context of stabilized term ð~ uc  r~ FEM can be found in [12,67]. With k the fluid iteration counter, the linearization of the convective term is given by [12,67]

~ v c;kþ1  r~ v kþ1  ~ v

c;k

 r~ v kþ1 þ k1~ v c;kþ1  r~ v k  k1~ v c;k  r~ v k:

ð32Þ

For k1 ¼ 0 the fixpoint and for k1 ¼ 1 the Newton–Raphson scheme is obtained. The present simulations are based on fixpoint iteration, because it possesses a more stable convergence rate. The Newton– Raphson scheme indeed possesses a higher convergence order, but

3326

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

it is quite sensitive to initial guesses. The flow simulations performed with the fixpoint scheme converged within 3–7 iterations. Thereby, no performance improvement is recognized with the Newton scheme. That is why the fixpoint scheme is utilized. 3.2. Computational solid mechanics The starting point of the spatial discretization with finite-elements is the weak (variational) formulation, which reads for (3) as follows: Find ~ u 2 ðH10 ðXs ÞÞd such that

~ Þ þ ðcB~ ~ ÞT Þ ¼ ð~ ~Þ qs ð~ u€; w u; ðBw f X; w ðH10 ð

ð33Þ

d

~2 for all w Xs ÞÞ . The spatial discretization then leads to the following semi-discrete Galerkin formulation

€ þ Ku ¼ F; Mu

ð34Þ

with M, K the finite-element matrices, u the nodal vector of unknown displacements, and F the right-hand side containing the fluid forces. In a similar way, the variational formulation for the geometrically nonlinear case can be obtained [4,41]. The application of the standard Galerkin FEM leads to a nonlinear matrix equation system

The integration by parts is also applied to the term including Lighthill’s tensor

! @ 2 T ij ;q @xi @xj

¼ Xf

    @T ij @T ij @q ni ; q  ; : @xj @xj @xi Xf Cfs

ð40Þ

Based on the continuity equation the term including @T ij =@xj can be substituted by

      @T ij @Pa @ qv i ;q ni ; q ¼  : ni ; q @xj @t @~ n Cfs Cfs Cfs

ð41Þ

The last term in (41) vanishes along an interface to a fixed solid. Along a flexible body, the following relation holds (see (22))

  @ qf v i ¼ ni ; q @t Cfs

@2u qf 2i ni ; q @t

! :

ð42Þ

Cfs

Using the relations (40) and (41) and combining the surface integrals of the acoustic pressure results in the variational formulation of flow- and vibration-induced sound

1 @ 2 Pa ;q c2 @t 2

! þ Xa

    @Pa @q 1 @Pa ; þ ;q c @t @xi @xi Xa Co |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}

Absorbing boundary

€ þ ðK þ Knlin ðuÞÞu ¼ F: Mu

ð35Þ

The semi-discrete Eqs. (34) and (35) are discretized in time with a Newmark scheme in effective stiffness form [38]. All presented simulations are based on Newmark integration parameters b ¼ 0:25 and c ¼ 0:5. The geometrically nonlinear system is finally linearized with a Newton–Raphson scheme [4].

!   @T ij @q @ 2 ui ¼ ;  qf 2 ni ; q : @xj @xi X @t Cfs |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl}f |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl}

ð43Þ

Flow-induced sound Vibrational-induced sound

The Galerkin FEM discretization finally yields the following DAE

€ a þ DP_ a þ KPa ¼ F; MP

ð44Þ

3.3. Computational acoustics Mathematically the fluid–acoustic interaction can be represented by Lighthill’s equation (17). In order to describe the finiteelement discretization of the fluid–acoustic coupling, we write this equation in index notation

3.4. Fluid–solid interaction

1 @ 2 Pa @ 2 Pa @ 2 T ij  ¼ ; 2 c2 @t2 @x @xi i @xj

ð36Þ

with appropriate initial and boundary conditions. The computational domain Xa is composed of a part Xf , on which the fluid mechanics quantities are of relevance, and a pure acoustic propagation part Xp ðXa ¼ Xf [ Xp Þ. At this stage, it is assumed that the flow variables are known inside the fluid domain Xf , e.g., due to a previously performed CFD simulation, and therefore the Lighthill tensor is known as well. To develop the variational form, (36) is multiplied by a test function q 2 H10 ðXa Þ and integrated over the whole domain Xa

! 1 @ 2 Pa @ 2 Pa @ 2 T ij   ; q ¼ 0: c2 @t 2 @xi @xj @x2i X

ð37Þ

a

The integral theorem of Green, applied to the spatial derivative of the acoustic pressure, yields

@ 2 Pa ;q @x2i

! ¼ Xa

    @P a @Pa @q  ; : ;q @xi @xi Xa @~ n Cfs [Co

ð38Þ

In (38) Cfs denotes the interface between fluid and solid, and Co the outer (limiting) boundary of the acoustic computational domain (see Fig. 2) at which we apply an absorbing boundary condition of first-order according to [41], expressed by

@Pa 1 @Pa ¼ : c @t @~ n

which is discretized in time with the Newmark scheme. As the impact of the acoustic quantity on the flow as well as on the vibrating structure is neglected, a sequential coupling scheme can be applied to resolve fluid–solid–acoustic coupling.

ð39Þ

As already mentioned, the fluid–solid interaction is a surface coupled problem and realized in a partitioned approach. In order to prescribe the coupling conditions, all elements need to be defined along the fluid–solid interface. The fluid forces according to (16) cause a deformation of the solid. The resulting solid field deformation leads to a domain update for the fluid field. In order to interpolate the deformation inside the fluid domain, a grid adaption is performed. The movement of the fluid domain is governed by the interface displacements ~ d. As soon as the interface displacements are larger than the fluid element size along the fluid–solid interface, it is necessary to interpolate the interface displacements inside the whole fluid domain in order to avoid overlapping elements. Hereby, the fluid elements should be deformed as little as possible to keep the numerical errors small. Different strategies exist to perform this grid adaption, e.g.,:  linear interpolation [26,56],  elliptic smoothing [70] and  the pseudo-structure approach [62,67]. The first two approaches can be implemented very efficiently. However, especially linear interpolation is restricted to simple deformations. Elliptic smoothing schemes already allow a more complex deformation state. The most powerful grid adaption scheme is given by the pseudo-structure scheme as it allows very complex deformations. Hereby, the fluid domain is assumed to be a

3327

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

pseudo-solid which is loaded by inhomogeneous Dirichlet boundary conditions at the interface, namely the interface displacements ~ d. The pseudo-domain deformation finally provides the grid adaption needed. If grid adaption fails, re-meshing may be necessary. However, the required projection of field variables from the old to the new mesh introduces projection errors and should therefore be avoided as long as possible [62]. The pseudo-solid approach is based on static solid mechanics where a pseudo-deformation field ~ r has to fulfill

with ~ r ¼~ d on Cfs :

r  rg ¼ 0 on Xf

ð45Þ

The pseudo stresses rg can be computed via the pseudo strains g ¼ 12 ðr~ r þ ðr~ rÞT Þ and the constitutive equation of a Hook body rg ¼ c g g . The domain of the grid mechanical field is generally the same as the fluid domain. As soon as the fluid domain deformation is bounded locally, it is possible to split the fluid domain Xf into a moving XALE and fixed part XEuler with

Xf ¼ XEuler [ XALE : In Fig. 3 such a local domain decomposition is sketched. Thereby, the grid adaption can be restricted to the ALE domain. The pseudo-solid field is discretized with the standard Galerkin FEM, yielding the LAS

Kr ¼ F:

ð46Þ

Here, the force vector F consists only of contributions from inhomogeneous Dirichlet boundary conditions at the fluid–solid interface Cfs . Furthermore, r denotes the vector of unknown grid displacements. A time discretization is hereby not needed. However, as the mesh velocities are required due to the ALE description of the fluid field, the grid velocities v g have to be computed based on the mesh displacements r. The following approximations of the grid acceleration ag and grid velocities v g are applied to obtain secondorder accurate time discretization

rnþ1  2rn þ rn1 ; Dt 2 nþ1 n r r Dt ¼  ag : 2 Dt

ag ¼

ð47Þ

vg

ð48Þ

In order to avoid badly shaped elements during grid adaption, the elasticity modulus varies from element to element. An exponential elasticity decrease provides good results for many applications. The elasticity modulus Edist of the pseudo-solid is thereby reversely proportional to the distance jx  xCfs j between a certain fluid element and the fluid–solid interface Cfs

 Edist ¼ E0

xmax j~ x ~ xCfs j

q :

ð49Þ

In (49) xmax denotes the distance between the fluid–solid interface and the fluid boundary. E0 stands for the initial elasticity modulus (assumed constant for all elements). Numerical tests showed that the best grid adaption is obtained with an exponent q ¼ 1:5 [48]. As shown in Fig. 4, the element elasticity Eelem is set reverse proportional to the distance between each element and the fluid–solid

Fig. 4. Grid adaption algorithm.

interface Cfs according to (49). In a second step, a strain-dependent contribution to the stiffness is computed, based on the pseudo strains g . The strain is computed and the elasticity is modified again in such a way that elements with a higher strain energy become even stiffer [10]

Estrain

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2g;11 þ 2g;22 : / 2

The pseudo-solid computation is now based on the composite elasticity distribution Eelem ¼ Edist þ Estrain . The stopping criterion is reached as soon as the maximal strain is below a prescribed bound. In order to speed-up or to make the fluid–solid convergence possible at all, the partitioned coupled scheme is treated with a Richardson iteration scheme [51]. Thereby, the following relaxation is introduced into the fluid–solid interaction nþ1 ^ nþ1 þ ð1  x Þdnþ1 ; dkþ1 ¼ xk d k k kþ1

ð50Þ

^ nþ1 . The relaxation with the un-relaxed interface displacement d kþ1 parameter x can thereby be chosen to be fixed or flexible during the simulation time. The flexible relaxation is according to Aitken [40] because it is known to be efficient and robust [51]. In each fluid–solid iteration k the relaxation parameter xk is computed based on the change of interface displacements Dd nþ1

nþ1

Ddkþ1 ¼ dk

^ nþ1 ; d kþ1

with the Aitken factor



ð51Þ

lnþ1 k

nþ1 lnþ1 ¼ lnþ1 k k1 þ lk1  1 



nþ1

Ddk

nþ1

 Ddkþ1 nþ1

D dk

T

nþ1

Ddkþ1

nþ1

 Ddkþ1

:

ð52Þ

The relaxation parameter computes thereby as

xnþ1 ¼ 1  lnþ1 k k :

ð53Þ

Alternative relaxation parameter estimations are, e.g., the reduced order model [64] and the gradient approach [51,67]. We want to note that instabilities were observed for our coupled fluid–solid scheme. The source of the instabilities is the socalled artificially added mass effect. In accordance to [23], we found the following relations:

Fig. 3. ALE and Euler domains of fluid.

 The instabilities occur earlier for decreasing time step sizes.  The instabilities get worse for a higher density ratio qf =qs as well as an increasing fluid viscosity.

3328

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

kdkþ1  dk k2 < efsi ; kdkþ1 k2

ð55Þ

whereby efsi is the defined tolerance. For further details we refer to [48]. 4. Validation examples The scheme presented has been validated by a number of test cases from which we will present two examples. The first one is a fluid–solid interaction and the second one an aeroacoustics test case. The coupling of solid and acoustics has already been validated by many computations [41]. 4.1. Fluid–solid interaction

Fig. 5. Fluid–solid–acoustic interaction scheme.

 The instabilities decrease for increasing solid stiffness.  The temporal discretization influences the onset of instabilities.

3.5. Fluid–solid–acoustic coupled algorithm Combining the fluid–solid interaction with the fluid–acoustic and solid–acoustic coupling finally yields the fluid–solid–acoustic scheme. The iterative coupling scheme to resolve fluid–solid interaction can be extended by the sequential coupling scheme of the fluid–acoustic and solid–acoustic couplings as shown in Fig. 5. To measure the convergence, the relative change of the interface displacements is tested by a stopping criterion. With ~ d as the solid displacements at the interface

~ d ¼~ u

on Cfs ;

the stopping criterion is given by

ð54Þ

To validate the scheme with respect to the fluid–solid coupling, we present computational results of the DFG benchmark setup displayed in Fig. 6. Concerning the experimental measurements, we refer to [28] and for simulations to [70]. As depicted in Fig. 6 a thin sheet (flag) is placed in the wake of a cylinder in a channel filled with a mixture of Polyglycol and water. The gravitational acceleration g is acting in x-direction. The maximal deflection of the trailing edge of the flag is measured to be approximately 20 mm in y-direction, so that the artificially added mass effect and mechanical geometric nonlinearities play a significant role. The clamping boundary conditions of the solid field are restricted to the origin of the cylinder, which is therefore able to rotate. The solid is composed of three materials and their material parameters are listed in Table 2. The fluid mechanical material parameters for the Polyglycol-water mixture are l ¼ 1:64  104 Pa s for the 3 dynamic viscosity and qf ¼ 1050 kg=m for the density. The time step size is set to Dt ¼ 5 ms for the fluid–solid simulation. As sketched in Fig. 6, a velocity-driven flow is considered with the following inflow profile

vx



jyj ¼ 1:45 1  120

30 !

m=s:

ð56Þ

At the top and the bottom boundary, wall conditions are applied and the outlet is realized with open boundary conditions. The stopping criterion of the incremental change of the fluid field is 1  103 m=s. The nonlinear solid is assumed to be converged if the error in the strain energy is less than 1  106 J. The stopping criterion of the absolute change of interface displacement of the fluid–solid interaction is set to 1  102 mm. The fluid–solid interaction simulation is restarted on a fluid simulation with a small y-velocity of 0.05 m/s at the flag in order to speed-up the onset of oscillations. The von Kármán vortex street in the wake of the

Fig. 6. Geometry of the benchmark problem for fluid–solid interaction (units in mm).

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334 Table 2 Material parameters. Material

E 11

Aluminum Steel 1 Steel 2

1:5  10 Pa 2:0  1011 Pa 2:0  1011 Pa

m

qs

0.3 0.3 0.3

4300 kg=m3 7855 kg=m3 7800 kg=m3

3329

of aeroacoustic codes ([19,52] and references therein). The inhomogeneous part of Lighthill’s equation is directly computed by the analytical formulas as provided by [52]. The configuration of the vortices is depicted in Fig. 8. The induced flow field is assumed to be inviscid and incompressible. The two point vortices are separated by a distance of 2r 0 and have a circulation Cv . The vortex cores rotate around each other with a period of T ¼ 8p2 r 20 =Cv . The angular speed of the swirling flow is given by x ¼ Cv =4pr 20 and the circumferential Mach number is Mar ¼ Cv =4pr 0 c0 . Since at the origin of the setup a singularity occurs, we apply a vortex core model in order to compute the flow velocity (for details, see [19]). The angular velocity of the vortex is set according to the circulation Cv ¼ 1:00531 m2 =s and the radius to r 0 ¼ 1 m. The speed of sound is set to c ¼ 1 m=s yielding a circumferential Mach number of Mar ¼ 0:08. The domain considered has an extension of 1000 m  1000 m in x- and y-direction. The numerical grid is sketched in Fig. 9a. In order to resolve the vortex source accurately, an orthogonal fine grid is clustered around the origin, with 300  300 Q2 elements. The adjoining unstructured mesh is

Fig. 7. DFG fluid–solid benchmark; ðx; yÞ-displacements of the trailing edge of the flag.

cylinder governs the flag deformation. The simulated and measured ðx; yÞ-displacements of the trailing edge correlate in an acceptable range, as shown in Fig. 7. It can be seen that the time step size for the simulation is quite large. Fluid–solid convergence problems are observed for smaller time step sizes which are traced back to the artificially added mass effect. The upper and lower bounds of the relaxation parameter is a crucial point for that problem. For the presented simulation, it is forced to be within the range of 0:01 6 xi 6 0:3. The differences between the simulation and the measurement are supposed to be based on the errors of the time discretization. 4.2. Aeroacoustics We consider the generation of acoustic sound by a rotating quadrupole, which is a commonly used test case for the validation

Fig. 8. Schematic drawing of co-rotating vortices.

Fig. 9. Fluid–acoustic test case: co-rotating vortex pair: (a) mesh (b) contour plots of the acoustic pressure in mPa distribution at t ¼ 400 s.

3330

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

Fig. 10. Decay of the acoustic pressure values along the x-axis.

composed of 108.000 second-order triangular elements which are embedded in a structured mesh of 250.000 Q2 elements. The fluid velocities are only prescribed analytically in the innermost rectangle ð6 m  6 mÞ. The Lighthill sources are taken into account. At the domain surrounded the sound propagation is computed with the homogeneous wave equation. Fig. 9b displays the contour plot of the computed acoustic pressure and Fig. 10 shows a quantitative comparison of the pressure distribution along the x-axis compared to the analytical solution. Overall, the acoustic field predicted is confirmed by the analytical solution. Due to the vortex core model, we have not displayed the numerical computation near the origin.

The mechanical field possesses clamped boundary conditions along the line where the vocal folds touch the trachea wall. Fluid–solid interaction conditions are used at the rest of the mechanical boundary. Linear elasticity is assumed for the mechanical field. Absorbing boundary conditions are prescribed at the inflow and outflow for the acoustic simulation. Inside the whole acoustic domain Lighthill’s source terms are evaluated according to the fluid–acoustic coupling. The time step size is set to Dt ¼ 0:1 ms. The 2D finite-element mesh consists of second-order quadrilateral elements. Details are shown in Fig. 12. The fluid and the acoustic domain are composed of 22.872 elements, the mechanical domain is discretized with 3.022 elements and the grid mechanical domain consist of 14.072 elements. The fluid field possess 226.074 DOFs (Degree-Of-Freedom), the solid 18.380 DOFs, the acoustic 69.441 DOFs and the grid adaption 84.014 DOFs. The fluid field converged within approximately 5 fluid iterations, and the fluid– solid interaction within 15 iterations in the mean. The computing time per deformation cycle is approximately 10 hours on four Intel Xeon 5160 processors of a Woodcrest computer-cluster (overall peak performance: 10.4 Tflops/s). The setup as reported in [63] was enhanced in our study to model the real physiology more accurately. Therefore, the model is divided into three sub-layers [2], the body, the ligament and the cover as shown in Fig. 13. The following abbreviations are used subsequently: vocal fold (vf), sub-glottal (sub), and supra-glottal (supra). The elasticity modulus of body, cover and ligament is 40,10 and 100 kPa. The Poisson number is set to 0.4 for each material. Under these elasticities the first vocal fold eigenfrequency is at 100 Hz. The simplification to the 2D plane setup has been affirmed with 3D simulations [48].

5. Fluid–solid–acoustic coupled simulations of human phonation 5.1. Phonation model The geometric assembly of the vocal folds is adopted to the model presented in [63]. The view of the larynx model in Fig. 11 is displayed rotated. Therefore, the air streams from the left to the right. The fluid domain is subdivided into two sub-domains: XALE and XEuler (see Fig. 11). In the domain XALE the fluid–solid interaction takes place. This part of the mesh has to be adapted within each fluid–solid iteration (ALE mesh). The rest of the fluid domain XEuler does not need any mesh adaption (Euler mesh).

Fig. 11. Computational larynx model (distances are given in mm).

Fig. 12. Mesh around the vocal fold: (a) zoom around the vocal folds; (b) zoom around the glottis.

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

3331

Fig. 13. Three layer vocal fold model (distances are given in mm).

5.2. Fluid–solid–acoustic coupled results In the following, two scenarios are considered in more detail. Initially, the development of the Coanda effect is explored based on a constant pressure drop along the glottis as it may arise within a human body due to the contraction of the lung. Subsequently, the acoustic impact of the Coanda effect is investigated. Thereby, it is of advantage to constrain the volume-induced sound source to a single frequency, which can be achieved by prescribing an oscillating pressure drop. The pressure drop is therefore set to DP ¼ 1:0 þ 0:2 sinð2pfp tÞ kPa with a frequency of 100 Hz [3].

Fig. 15. Development of the Coanda effect within three cycles.

5.2.1. Development of the Coanda effect In order to investigate the development of the Coanda effect, a constant pressure drop of 500 Pa is adjusted between the inlet and the outlet. The initial glottis width is set to dGlottis ¼ 0:5 mm. The mechanical deformations of the fluid-induced vibrations for one deformation cycle is displayed in Fig. 14. Starting from state (a) the vocal folds deform and tend to close the gap and the glottis becomes convergent (b). However, at the turning point (c) the glottis shape becomes divergent (d). During the back movement of the vocal folds the shape again becomes convergent (e to f). This deformation cycle is repeated continuously. The presented dynamics coincide with experiments using excised human larynges [14,15].

Fig. 16. Mechanical displacements in streamwise direction (vertical/x-direction) of point A (Fig. 11) for a time span of 0.03 s.

Fig. 15 captures the time span of the first four cycles with contour plots of the velocity magnitude around the deformed vocal folds. The flow boundary conditions are basically symmetric and therefore a symmetric jet develops first. Owing to the instability of the developing shear layer (Kelvin–Helmholtz instability) jet oscillations appear downstream with increasing simulation time. This instability moves upstream in time and the corresponding oscillations further increase. Finally, the jet attaches to the trachea wall. The side to which the jet attaches is thereby randomly.

Fig. 14. Computed deformation cycle of the vocal folds (vf).

5.2.2. Acoustic impact of the Coanda effect The main focus is now to investigate the impact of the Coanda effect on the sound generation. Between the inflow and outflow a oscillating pressure drop as described above is assumed to constrain the volume-induced sound source to a single frequency. The frequency of the pressure drop oscillation is chosen to be

3332

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

Fig. 19. Fluid velocity at point B (Fig. 11) in streamwise direction (vertical/xdirection) of the asymmetric case, i.e. with Coanda effect.

Fig. 17. Velocity and nodal acoustic loads: (a) velocity magnitude with Coanda effect; (b) velocity magnitude under symmetric flow conditions; (c) normalized nodal acoustic loads due to (a); (d) normalized nodal acoustic loads due to (b).

Fig. 20. Fluid velocity at point D (Fig. 11) in cross-stream direction (lateral/ydirection) of the asymmetric case, i.e. with Coanda effect.

Fig. 18. Mechanical displacements at point A (Fig. 11) in streamwise direction (vertical/x-direction) of the asymmetric case, i.e. with Coanda effect.

fp ¼ 100 Hz. The cross-stream velocity at the inlet is set to zero and open boundary conditions are defined at the outlet. Wall no-slip/ no-penetration boundary conditions are prescribed at the trachea and at the interface to the vocal folds fluid–solid interaction.

The history of the streamwise displacement of point A (see Fig. 11) is shown in Fig. 16 for a time slot including three deformation cycles. These displacement values correspond to earlier experimental results [6]. The mechanical deformation thereby oscillates with 100 Hz. The complete simulation time is set to 0.36 s comprising 36 cycles and a frequency resolution of 2.8 Hz is obtained, when a Fourier transformation is applied. A typical velocity magnitude distribution is shown in Fig. 17a. The air jet downstream of the glottis is highly unsteady. The wall varies thereby, at which the air jet attaches. The finite-element discretization of Lighthill’s acoustic analogy provides nodal acoustic loads, representing acoustic point sources [19]. In Fig. 17c the spatial distribution of the nodal acoustic loads due to the velocity field of Fig. 17a is displayed. Especially in the shear layers, large values of the acoustic nodal loads are observed as expected. In a second computation, we enforce symmetric flow conditions throughout the whole simulation time and arrive at a hemilarynx model, which will show no Coanda effect. The geometry is the same as shown in Fig. 11, however the hemilarynx model considers only one half of the domain. At the midline symmetry conditions are applied. Under symmetric flow conditions a typical velocity distribution and the corresponding nodal acoustic loads are displayed in Fig. 17b/d. The main sound generating region of the

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

3333

on the air jet fluctuations, which are highly unsteady. The influence of the jet fluctuations on the acoustic pressure can be seen in Fig. 21a. The 100 Hz peak is again dominant, but strong lower and higher frequency contributions can be recognized. It can be concluded that the acoustic pressure is based on the complete spatially distributed air jet fluctuations. Under symmetric flow conditions of the hemilarynx model, there is no cross-stream velocity at point D and the acoustic pressure spectrum (Fig. 21b) contains just the dominant 100 Hz peak without other significant frequency components. The frequency spectra of the streamwise displacement and the streamwise velocity of the symmetric case are omitted, because these quantities have the same frequency spreading as in the asymmetric case (Figs. 18 and 19). Remarkable is that the base signal is strongly enriched especially with lower frequency components, when the Coanda effect and thus a strongly asymmetric flow field can develop. The Coanda effect during phonation has also be found in former studies [60,18,35]. However, its influence on voice quality is still not fully understood.

6. Conclusions The objective of this contribution was to develop a computational scheme resolving the fluid–solid–acoustic coupling. Many investigations considering either the fluid–solid interaction or the fluid–acoustic interaction have been undertaken. These two scientific branches were coupled within this contribution. The application of this combined scheme to investigate human phonation has been demonstrated and provided new results, especially to the impact of the Coanda effect and the resulting highly asymmetric and instanteneous flow field on the frequency spectra of the acoustic sound. Furthermore, the presented simulations showed realistic self-sustained vocal fold vibrations with converging shape during the opening and diverging shape during the closing phase as in-vivo dynamics. Currently, we extend our scheme to visco-elastic material modeling of the vocal folds. Furthermore, we plan to include mechanical contact for the vocal folds. Fig. 21. Power spectrum of the acoustic pressure at point C (Fig. 11) (a) with a developed Coanda effect; (b) under symmetric flow conditions, i.e. without Coanda effect.

symmetric flow conditions is near the jet front and the shear layers, while by asymmetric conditions the region where the jet attaches the wall has a significant contribution to the sound. Furthermore, the region just downstream the glottis is a noteworthy sound source (Fig. 17c). In both regions strong shear layers exist. The nodal acoustic loads, and therefore also the acoustic field, depend strongly on the air jet dynamics for a developed Coanda effect. This air jet dependency can be exhibited by frequency spectra of several field variables. At first, the frequency spectra of the asymmetric case are discussed followed by the symmetric case. In the frequency spectra of the mechanical displacements at point A (Fig. 18) and in the fluid velocity at point B (Fig. 19), both in streamwise direction, the 100 Hz peak is dominant and no other significant frequency components are contained. Both the streamwise displacement and the streamwise velocity are mainly dependent on the prescribed pressure oscillation. However, the cross-stream velocity of point D (see Fig. 11) gives an impression of the air jet fluctuations. In the frequency spectrum of the cross-stream velocity of point D (Fig. 20) the 100 Hz component does not exist. The cross-stream velocity is mainly dependent

Acknowledgements This work was supported by Deutsche Forschungsgemeinschaft Grant No. FOR894 Strömungsphysikalische Grundlagen der menschlichen Stimmgebung. References [1] R.A. Adams, Sobolev Spaces, Academic Press, 1975. [2] F. Alipour, D.A. Berry, I.R. Titze, A finite-element model of vocal-fold vibration, J. Acoust. Soc. Am. 108 (2000) 3003–3012. [3] F. Alipour, R.C. Scherer, Flow separation in a computational oscillating vocal fold model, J. Acoust. Soc. Am. 116 (2004) 1710–1719. [4] K.J. Bathe, Finite-element Procedures, Prentice Hall, 1996. [5] C. Bogey, C. Bailly, D. Juvé, Computation of flow noise using source terms in linearized Euler’s equations, AIAA J. 40 (2) (2002) 1235–1243. [6] A. Bößenecker, D.A. Berry, J. Lohscheller, U. Eysholdt, M. Dllinger, Mucosal wave properties of a human vocal fold, Acta Acust. United Acust. 93 (5) (2007) 815–823. [7] F. Brezzi, M.O. Bristeau, L.P. Franca, M. Mallet, G. Rogé, A relationship between stabilized finite-element methods and the Galerkin method with bubble functions, Comput. Meth. Appl. Mech. Engrg. 96 (1992) 117–129. [8] F. Brezzi, M. Fortin, Mixed and Hybrid Finite Element Methods, Springer, 1991. [9] A.N. Brooks, T.J.R. Hughes, Streamline upwind/Petrov–Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Engrg. 32 (1982) 199–259. [10] G. Chiandussi, G. Bugeda, E. Onate, A simple method for automatic update of finite-element meshes, Commun. Numer. Meth. Engrg. 16 (2000) 1–19. [11] H. Coanda, Device for deflecting a stream of elastic fluid projected into an elastic fluid US-patent, 2.052.869, 1935.

3334

G. Link et al. / Comput. Methods Appl. Mech. Engrg. 198 (2009) 3321–3334

[12] R. Codina, A nodal-based implementation of a stabilized finite element method for incompressible flow problems, Int. J. Numer. Meth. Fluids 33 (2000) 737– 766. [13] D.G. Crighton, Computational aeroacoustics for low Mach number flows, in: J.C. Hardin, M.Y. Hussaini (Eds.), Computational Aeroacoustics, SpringerVerlag, 1992. [14] M. Döllinger, D.A. Berry, Visualization and quantification of the medial surface dynamics of an excised human vocal fold during phonation, J. Voice 20 (3) (2006) 401–413. [15] M. Döllinger, F. Rosanowski, U. Eysholdt, J. Lohscheller, Grundlagenuntersuchungen für Stimmlippendynamik: 3DSchwingungsanalyse menschlicher und caniner Kehlköpfe, HNO 56 (12) (2008) 1213–1220. [16] J. Donea, A. Huerta, Finite-element Methods for Flow Problems, Wiley, 2003. [17] B.D. Erath, M.W. Plesniak, The occurence of the Coanda effect in pulsatile flow through static models of the human vocal folds, J. Acoust. Soc. Am. 120 (2) (2006) 1000–1011. [18] B.D. Erath, M.W. Plesniak, The occurrence of the Coanda effect in pulsatile flow through static models of the human vocal folds, J. Acoust. Soc. Am. 120 (2006) 1000–1011. [19] M. Escobar, Finite-element Simulation of Flow-induced Noise Using Lighthill’s Acoustic Analogy, Ph.D Thesis, University of Erlangen-Nuremberg, 2007. [20] R. Ewert, W. Schröder, Acoustic perturbation equations based on flow decomposition via source filtering, J. Comput. Phys. 188 (2) (2003) 365–398. [21] J.H. Ferziger, M. Peric´, Computational Methods for Fluid Dynamics, Springer, 2002. [22] J.E. Ffowcs-Williams, D.L. Hawkings, Sound radiation by turbulence and surfaces in arbitrary motions, Philos. Trans. Roy. Soc. London 264 A (1969) 321–342. [23] Ch. Förster, W.A. Wall, E. Ramm, Artificial added mass instabilities in sequential staggered coupling of nonlinear structures and incompressible viscous flows, Comput. Meth. Appl. Mech. Engrg. 196 (7) (2007) 1278–1293. [24] L.P. Franca, C. Farhat, Bubble functions prompt unusual stabilized finiteelement methods, Comput. Meth. Appl. Mech. Engrg. 123 (1995) 299–308. [25] L.P. Franca, S.L. Frey, Stabilized finite-element methods: II The incompressible Navier–Stokes equations, Comput. Meth. Appl. Mech. Engrg. 99 (1992) 209– 233. [26] M. Glück, Ein Beitrag zur numerischen Simulation von Fluid-StrukturInteraktionen – Grundlagenuntersuchungen und Anwendung auf Membrantragwerke (A contribution to numerical simulation of fluidstructure interactions – Fundamental study and application to membranous structures), Ph.D Thesis, University of Erlangen-Nuremberg, 2003. [27] M.E. Goldstein, Aeroacoustics, McGraw-Hill, New York, 1976. [28] J.P. Gomes, H. Lienhart, Experimental study on a fluid-structure interaction reference test case, in: Fluid-Structure Interaction – Modelling, Simulation, Optimization, Springer, 2006. [29] V. Gravemeier, The Variational Multiscale Method for Laminar and Turbulent Incompressible Flow, Ph.D Thesis, University of Stuttgart, 2003. [30] P.M. Gresho, R.L. Sani, Incompressible Flow and the Finite-element Method, vol. 2, Wiley, 2000. [31] H.E. Gunter, A mechanical model of vocal-fold collision with high spatial and temporal resolution, J. Acoust. Soc. Am. 113 (2003) 994–1000. [32] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II - Stiff and Differential-Algebraic Problems, Springer Series in Computational Mathematics, second rev. ed., vol. 14, Springer, 2002. [33] W. Heller, J. Klingenberg. Beiträge zur Strömungsmechanik (Contributions to fluid mechanics), Technical University of Dresden, 2001. [34] C. Hirt, A.A. Amsden, J. Cook, An arbitrary Lagrangian–Eulerian computing method for all flow speeds, J. Comput. Phys. 14 (1974) 227–253. [35] G.C.J. Hofmans, G. Groot, M. Ranucci, G. Graziani, A. Hirschberg, Unsteady flow through in-vitro models of the glottis, J. Acoust. Soc. Am. 113 (2003) 1658– 1675. [36] T.J.R. Hughes, Multiscale phenomena: Green’s functions, the Dirichlet-toNeumann formulation, subgrid scale models, bubbles and the origins of stabilized methods, Comput. Meth. Appl. Mech. Engrg. 127 (1995) 387–401. [37] T.J.R. Hughes, W.K. Liu, T. Zimmermann, Lagrangian–Eulerian finite-element formulation for compressible viscous flows, Comput. Meth. Appl. Mech. Engrg. 29 (1981) 329–349. [38] T.J.R. Hughes, The Finite-Element Method, Prentice-Hall, New Jersey, 2000. [39] T.J.R. Hughes, G. Scovazzi, L.P. Franca, Encyclopedia of Computational Mechanics: Ch 2 Multiscale and Stabilized Methods, vol. 3, John Wiley & Sons Ltd., 2004. [40] B.M. Irons, R.C. Tuck, A version of the Aitken accelerator for computer iteration, Int. J. Numer. Meth. Engrg. 1 (1969) 275–277. [41] M. Kaltenbacher, Numerical Simulation of Mechatronic Sensors and Actuators, second ed., Springer, 2007. [42] M. Kaltenbacher, M. Escobar, S. Becker, I. Ali, Computational aeroacoustics based on Lighthill’s analogy, in: S. Marburg, B. Nolte (Eds.), Computational Acoustics of Noise Propagation in Fluids, Springer, 2008. [43] M.H. Krane, Aeroacoustic production of low-frequency unvoiced speech sounds, J. Acoust. Soc. Am. 118 (2005) 410–427.

[44] L. Lee, R.J. LeVeque, An immersed interface method for incompressible Navier– Stokes equations, SIAM J. Sci. Comput. 25 (3) (2003) 832–856. [45] S. Lele, Computational aeroacoustics: a review, AIAA J. (1997) 18. [46] M.J. Lighthill, On sound generated aerodynamically. Part I: general theory, Proc. Roy. Soc. London A211 (1952) 564–587. [47] M.J. Lighthill, On sound generated aerodynamically. Part II: turbulence as a source of sound, Proc. Roy. Soc. London A222 (1954) 1–32. [48] G. Link, A Finite-Element Scheme for Fluid–Solid–Acoustics Interactions and its Application to Human Phonation. Ph.D Thesis, University of ErlangenNuremberg, 2008. [49] H. Luo, R. Mittal, X. Zheng, S.A. Bielamowicz, R.J. Walsh, J.K. Hahn, An immersed-boundary method for flow-structure interaction in biological systems with application to phonation, J. Comput. Phys. 227 (2008) 9303– 9332. [50] R. Mittal, G. Iaccarino, An immersed boundary methods, Ann. Rev. Fluid Mech. 37 (1) (2005) 239–261. [51] D.P. Mok, Partitionierte Lösungsansätze in der Strukturdynamik und der FluidStruktur-Interaktion (Partitioned analysis schemes in structural dynamics and fluid-structure-interaction). Ph.D Thesis, University of Stuttgart, 2001. [52] E. Müller, F. Obermeier, The spinning vortices as a source of sound, CP-22 AGARD (1967) 22.1–22.8. [53] A. Oberai, F. Ronaldkin, T. Hughes, Computation of trailing-edge noise due to turbulent flow over an airfoil, AIAA Aeroacoust. J. 40 (2002) 2206–2216. [54] C.S. Peskin, An immersed boundary method, Acta Numer. 11 (1) (2002) 479– 517. [55] M.O. Rosa, J.C. Pereira, M. Grellet, A. Alwan, A contribution to simulating a three-dimensional larynx model using the finite-element method, J. Acoust. Soc. Am. 114 (5) (2003) 2893–2905. [56] F. Schäfer, S. Kniesburges, T. Uffinger, S. Becker, J. Grabinger, G. Link, M. Kaltenbacher, Numerical simulation of fluid-structure- and fluid-structureacoustic interaction based on a partitioned coupling scheme, in: High Performance Computing in Science and Engineering Munich 2007, Transactions of the Third Joint HLRB and KONWIHR Result and Reviewing Workshop, Leibniz-Rechenzentrum München (LRZ), Garching, Springer, 2007. [57] O. Schenk, K. Gärtner, Solving unsymmetric sparse systems of linear equations with PARDISO, J. Future Generat. Comput. Syst. 20 (3) (2004) 475–487. [58] J. Seo, Y. Moon, Linearized perturbed compressible equations for low Mach number aeroacoustics, J. Comput. Phys. 218 (2) (2006) 702–719. [59] C. Tao, J.J. Jiang, Y. Zhang, Simulation of vocal fold impact pressures with a selfoscillating finite-element model, J. Acoust. Soc. Am. 119 (6) (2003) 3987–3994. [60] C. Tao, Y. Zhang, D.G. Hottinger, J.J. Jiang, Asymmetric airflow and vibration induced by the Coanda effect in a symmetric model of the vocal folds, J. Acoust. Soc. Am. 122 (2007) 2270–2278. [61] T.E. Tezduar, S. Mittal, S.E. Ray, R. Shih, Incompressible flow computations with stabilized bilinear and linear equal-order-interpolation velocity–pressure elements, Comput. Meth. Appl. Mech. Engrg. 95 (1992) 221–242. [62] T.E. Tezduyar, Encyclopedia of Computational Mechanics: Ch. 17 Finiteelement Methods for Fluid Dynamic with Moving Boundaries and Interfaces, John Wiley & Sons Ltd., 2004. [63] S.L. Thomson, L. Mongeau, S.H. Frankel, Physical and numerical flow-excited vocal fold models, in: Third International Workshop MAVEBA, Firenze University Press, 2003, ISBN 88-8453-154-3, pp. 147–150. [64] J. Vierendeels, K. Dumontb, P.R. Verdonckb, A partitioned strongly coupled fluid-structure interaction method to model heart valve dynamics, J. Comput. Appl. Math. 218 (2008) 602–609. [65] M.P. de Vries, H.K. Schutte, A.E.P. Veldman, G.J. Verkerke, Glottal flow through a two-mass model: Comparison of Navier–Stokes solutions with simplified models, J. Acoust. Soc. Am. 111 (4) (2002) 1847–1853. [66] C. Wagner, T. Hüttl, P. Sagaut (Eds.), Large-Eddy Simulation for Acoustics, Cambridge University Press, New York, 2007. [67] W.A. Wall. Fluid-Struktur-Interaktion mit stabilisierten Finiten Elementen (Fluid-structure interaction with stabilized finite-elements). Ph.D Thesis, University of Stuttgart, 1999. [68] W.A. Wall, A. Gerstenberger, P. Gammnitzer, C. Förster, E. Ramm, Large deformation fluid-structure interaction – advances in ALE methods and new fixed grid approaches, in: H.-J. Bungrartz, M. Schäfer (Eds.), Fluid-Structure Interaction: Modeling, Simulation, Optimization, LNCSE, Springer, 2006. [69] A. Gerstenberger, W.A. Wall, An extended finite-element method/Lagrange multiplier based approach for fluid-structure interaction, Comput. Meth. Appl. Mech. Engrg. 197 (2008) 1699–1714. [70] S. Yigit, M. Schäfer, M. Heck, Grid movement techniques and their influence on laminar fluid structure interaction computations, J. Fluids Struct. 24 (6) (2008) 14. [71] W. Zhao, S.H. Frankel, L. Mongeau, Numerical simulation of sound from confined pulsating axisymmetric jets, AIAA J. 39 (10) (2001) 1868–1874. [72] W. Zhao, C. Zhang, S.H. Frankel, L. Mongeau, Computational aeroacoustics of phonation, Part I: computational methods and sound generation mechanisms, J. Acoust. Soc. Am. 112 (2002) 2134–2146. [73] Z. Zhang, L. Mongeau, S.H. Franke, S.L. Thomson, Sound generation by steady flow through glottis-shaped orifices, J. Acoust. Soc. Am. 116 (2004) 1720– 1728.