Contents lists available at ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

Numerical simulation of glottal ﬂow A. Hundertmark-Zaušková a,n, R. Lehmann b, M. Hess c, F. Müller c a b c

Institute of Mathematics, Johannes-Gutenberg University, Mainz, Germany Institute of Geosciences, Johannes-Gutenberg University, Mainz, Germany Department of Voice, Speech and Hearing Disorders, University Medical Center Hamburg-Eppendorf, Germany

art ic l e i nf o

a b s t r a c t

Article history: Received 4 October 2012 Accepted 10 September 2013

In cases of permanent immobility of both vocal folds patients have difﬁculties with breathing but rarely with voicing. However, clinical experience shows that the shape of the larynx (voice box) seems to have a signiﬁcant inﬂuence on the degree of airﬂow and breathing pattern. In order to ﬁnd an optimal geometry of the larynx in terms of easiness for breathing after the surgical change of vocal folds or false vocal cords (ventricular folds), a set of numerical simulations of glottal ﬂow for weakly compressible Navier–Stokes equations has been performed. We compare airﬂow resistance and volumetric ﬂow rate for several geometry concepts for inspiration as well as expiration. Finally, we discuss the optimal geometry with respect to the quality of breathing. & 2013 Elsevier Ltd. All rights reserved.

Keywords: Glottal ﬂow Numerical simulation Weakly compressible ﬂow Geometry optimization Respiration Airﬂow resistance Volumetric ﬂow rate

1. Respiration process and speaking In the respiration process (inspiration and expiration) and during speaking, air ﬂows from the lungs downstream through the windpipe (trachea), passes the vocal folds and the false vocal folds (ventricular folds) and leaves the larynx at the epiglottic cartilage into the pharynx. During respiration the vocal folds are actively held in a separated conﬁguration by muscular action. The elastic, vibrating vocal folds are responsible for voice production (so-called phonation). Position and size of the vocal folds differ between respiration (vocal folds are separated) and voice production (vocal folds are very close to each other but allow airﬂow). For voice production (i.e., phonation) the vocal folds are brought together actively by muscular actions. This conﬁguration of adducted vocal folds enables air pressure to build up under the closed (adducted) vocal folds until they are pushed apart due to their ﬂoppy, elastic tissue property. The aerodynamic-viscoelastic properties induce temporary closure of the vocal folds followed by immediate opening because of new air pressure build-up. An oscillatory cyclic behavior is started and makes the vocal folds vibrate regularly. The frequency of vocal fold oscillation is about 100 times per second in adult males, and approximately 200 times in females. The periodic pressure changes in the area of vocal folds generate sound. In this process air is ﬂowing gustily through the

n

Corresponding author. Tel.: þ 49 6131 3925099. E-mail address: [email protected] (A. Hundertmark-Zaušková).

0010-4825/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compbiomed.2013.09.009

vocal folds gap, which is about 1 mm wide, whereas it opens to about 10 mm during the inspiratory process, see Figs. 1 and 2. The laryngeal aerodynamics, the voice production (phonation) and their interaction have been extensively studied in laboratory experiments as well as using computational simulations in recent years. From the number of works we refer to, e.g., Scherer et al. [1–3,16,21,22], Khosla and Scherer [14], Bailly et al. [5], Apostoli [4], Titze et al. [23,24] and the citations therein. In these works ﬂow regimes, ﬂow vorticity, pressure drops and the impact of the surrounding larynx geometry on the voice generation have been studied. The impact of the ventricular constriction due to the false vocal folds has been clearly conﬁrmed. More precisely, the ventricular bands produce an additional pressure drop and affect the vocal folds closure which improves the voice quality [4,5]. In particular in pathological cases of human phonation or throat-singing phonation the mechanics of slower oscillations of false vocal folds play an important role. Furthermore, the shape of the vocal folds also signiﬁcantly affects the pressure distributions on the glottis and needs to be well speciﬁed when building computational and physical models, see [16,21,22]. For other works on modeling of the phonation process (with oscillating vocal folds) see, e.g., [17,19,20]. In [17,20] the ﬂuid– structure interaction problem is solved. The ﬂow development behind the non-vibrating rigid vocal folds has been studied for example in [6] for an incompressible and in [18] for a compressible viscous ﬂuid. Langnickel [15] belongs to the ﬁrst studies about the inﬂuence of the changed vocal folds’ shape on the airﬂow and its resistance. In Fig. 3 a possible conﬁguration of the gap between the vocal folds is shown in two dimensional coronal view. However, under

2178

A. Hundertmark-Zaušková et al. / Computers in Biology and Medicine 43 (2013) 2177–2185

Fig. 1. Laryngoscopic view of the vocal folds: open position for respiration (left), closed position for speaking (right), Medical Healthcom Ltd., Prague.

Fig. 3. Post-surgical situation with modiﬁed distances T, D>, A.

Fig. 2. The coronal cut of the larynx, Gray, Anatomy of the Human Body [9].

physical effort some patients complain about breathing difﬁculties (shortness of breath) causing the feeling of lack of air in the lungs. If vocal folds do not actively move (separation or adduction), the ‘ﬁxed’ gap between both vocal folds determines how breathing and voicing is. The more they are separated, the better the condition for respiratory airﬂow. Breathing is easier but voicing reduced. On the other hand, a very narrow gap results in good voicing but very limited breathing ability.

The aim of this work is to simulate the glottal ﬂow during the respiration process in vitro and to design a suitable geometry with respect to the air volume and ﬂow resistance after passing the modiﬁed vocal folds. Here, we assume the left vocal fold to be ﬁxed to the speaking position. This ensures the ability to generate sound as the right vocal fold can still elastically deform. However, the enlargement of the left vocal fold affects the breathing process, see [15]. Starting from this, our goal is to improve the breathing process by changing the shape of the larynx. For this purpose we design a simpliﬁed geometry of the larynx, see Fig. 3, and perform a set of experiments varying the distances T, D, A: T—distance of false vocal folds, D—thickness of the false vocal fold, A—distance between the vocal fold and the false vocal fold. Finally, a suitable geometry is proposed after testing several geometry conﬁgurations. We would like to point out that no optimization techniques have been applied. In this study the left non-vibrating

A. Hundertmark-Zaušková et al. / Computers in Biology and Medicine 43 (2013) 2177–2185

2179

vocal fold is enlarged (and ﬁxed to the speaking position), the right vocal fold is in opened position. The voice production based on models for acoustic waves and fold oscillation as well as the geometry optimization during the phonation is beyond the topic of this paper. This may be the goal of our future research.

310 K) pa ¼ 101325 Pa and air density ρa ¼ 1:138 kg=m3 , hence c 8:5 104 Pa=ðkg m 3 Þγ . Consequently, the density can be expressed using the pressure as follows: p þ ðpa pref Þ 1=γ ρ¼ : ð2:2Þ c

2. Mathematical model and computational data

Note that p denotes the computed pressure deviation from the atmospheric pressure. Since we start our simulation with an initial reference pressure of pref ¼ 100 Pa, we added pa pref to p in the above relation for the realistic representation of the density. Now the solution of (2.1) reduces to the pressure (deviation) p and the velocity v.

2.1. Governing equations We consider a simpliﬁed two dimensional coronal cut of the larynx, see Figs. 2–4. The vocal folds are assumed to be rigid. Since the occurring (air) velocities inside the larynx are much smaller than the speed of sound, the Mach number is small and, therefore, the air ﬂow can be modeled using the weakly compressible Navier– Stokes equations. The additional assumption of adiabatic and barotropic ﬂow leads to the independency of the conservation of energy from the system; thus, the conservation of energy is not considered in our model. According to the barotropic relation, the density of the ﬂuid ρ is given by the pressure p via p ¼ cργ , where c is a constant (we address to it later) and γ is the isentropic exponent, γ ¼ 1:4 for air. The resulting equations read as follows:

ρt þdivðρvÞ ¼ 0;

ð2:1aÞ

ðρvÞt þ divðρv vÞ ¼ ρf ∇p þ divð2μDðvÞ λ div vÞ

ð2:1bÞ

p ¼ cργ :

ð2:1cÞ

Here v denotes the velocity, ρ density, p pressure, μ; λ viscosity coefﬁcients and DðvÞ ¼ 12ð∇v þ ∇vT Þ the symmetric stress tensor. In our simulations we use the following values. The dynamic viscosity μ ¼ 1:8 10 5 Pa s, λ ¼ 23μ, [8,13]. The constant c in (2.1c) can be determined from the atmospheric pressure (with temperature

2.2. Boundary and initial conditions We use the following initial conditions for pressure and velocity in the whole computational domain Ω, see Fig. 4: v0 ðxÞ ¼ vðx; 0Þ ¼ 0 m=s;

xAΩ

p0 ðxÞ ¼ pðx; 0Þ ¼ pref ¼ 100 Pa;

x A Ω:

On the walls of the trachea Γ , see Fig. 4, we prescribe the non-slip boundary condition: vðx; tÞ ¼ 0

for x A Γ ;

t 4 0:

For the inspiration process we used the following boundary conditions at the remaining boundaries. At the top boundary Γ top (larynx– pharynx passage) we consider Dirichlet boundary conditions for the velocity; the horizontal velocity is set to zero, the vertical component is constant with respect to the spatial coordinate x and given by a time-dependent sine with an amplitude of 4 m/s: v1 ðx; tÞ ¼ 0;

π v2 ðx; tÞ ¼ 4 sin t 2

for x A Γ top ;

Fig. 4. The simpliﬁed geometry of the coronal cut of the larynx, breathing (left), speaking (right).

0 o t o 2 s:

2180

A. Hundertmark-Zaušková et al. / Computers in Biology and Medicine 43 (2013) 2177–2185

Fig. 5. Inspiration (t¼ 1 s), left: the velocity (m/s), right: the pressure distribution (Pa).

The amplitude 4 m=s corresponds to the maximum inspiration velocity at the time t ¼ 1 s, at the entrance to the human larynx, approximately ﬁtting the tidal volume of one breath at moderate exercise (which is 2.2 l for a man), see [12].1 At the bottom boundary Γ bottom (trachea) we consider the so-called ‘do nothing’ condition, i.e., zero normal stress during inspiration: ð pI þ μDðvÞ 23μ div vIÞ n ¼ 0

for x A Γ bottom ;

0 o t o 2 s:

For the expiration process the boundary conditions are swapped. At the bottom boundary Γ bottom we consider the same Dirichlet boundary condition for velocity as for Γ top above with opposite sign: π for x A Γ bottom ; 0 o t o2 s: v1 ðx; tÞ ¼ 0; v2 ðx; tÞ ¼ 4 sin t 2 On the outﬂow boundary which is now Γ top we ﬁx the pressure to the initial value of 100 Pa and zero the viscous stress μDðvÞ 23μ div v I ¼ 0. For our simulations we used the software COMSOL Multiphysics. The implemented ﬁnite element method allows to use several types of elements (Lagrange polynomial of the ﬁrst, second and third orders), for instance the P1–P1, P2–P1 or P3–P2 elements. The presented results are obtained using ﬁrst order Lagrange polynomial for pressure as well as velocity (P1–P1). For time discretization the ﬁrst and the second order backward differentiation formula is used. For more details see [7].

3. Experiments and measurements We start this section with simulations in a healthy human larynx presented by the simpliﬁed two dimensional cut, see Fig. 4 (left). Figs. 5 and 6 show the streamlines, the velocity ﬁeld and the pressure during breathing at the time instance t ¼ 1 s. This time instance represents the middle of the inspiration and expiration cycles, when the velocity achieves its maximum.

We introduce some physical quantities for measuring the quality of breathing for later comparison of the original status and the modiﬁed larynx geometries. In what follows we intend to ﬁnd the optimal geometry with respect to these quantities. In order to measure the volume of air ﬂowing into the lungs, R we use the ﬂow rate Q~ ¼ Γ v nΓ ; where Γ ¼ Γ bottom is the bottom boundary, see Fig. 4, and nΓ is the outward-pointing unit normal to this boundary. Note that nΓ ¼ ð0; 1ÞT . As we consider the longitudinal section of the trachea, the boundary Γ bottom is a one dimensional line. Aiming to compute realistically dimensionalized values, we approximate the elliptic cross-sectional shape of the trachea using a rectangle of the same ~ ~ area with edges a, b. We assume a ¼ 1:25 cm and b ¼ 0:7a ¼ 0:875 cm as elliptic semiaxes, see, e.g., [9], yielding an area of 3:44 cm2 . The length of the outﬂow boundary Γ bottom equals the transverse elliptic diameter a~ ¼ 2a ¼ 2:5 cm, thus, we obtain an approximation of the volumetric ﬂow rate Q: Z v nΓ bottom : Q b~ Q~ ¼ b~ Γ bottom

The ﬂow rate is computed in liters per second. During the expiration it is useful to evaluate the ﬂow rate above the vocal cords, i.e., the developed ﬂow on the top boundary Γ top . Note that the Dirichlet boundary condition for velocity at the bottom strongly affects the ﬂow rate there. Moreover, we compute the pressure difference Δp between Γ top , Γ bottom in the larynx domain with the help of mean pressures: 1 Z Δp ¼ jptop pbottom j ¼ pðx; tÞ dx jΓ top j Γ top

pðx; tÞ dx; jΓ bottom j Γ bottom 1

Z

t ¼ 1 s;

where the lengths jΓ top j and jΓ bottom j of the boundaries are 3 cm and 2.5 cm, respectively. We deﬁne the airﬂow resistance R ¼ Δp=Q

1

Values 2.2 l (male) and 1.4 l (female) are given in [12, Table 3]. The boundary velocity 4 sin ðπ t=2Þ m=s integrated over ½0; 2 s and multiplied with elliptic area of larynx–throat passage with semiaxes 1.5 cm, 1.05 cm yields 2:5 l. Thus, using the amplitude of 4 m/s gives an approximation to the male tidal breath volume.

relating the pressure difference to the ﬂow rate, cf. [10,15]. For optimization we will focus on the ﬂow rate Q and the airﬂow resistance R.

A. Hundertmark-Zaušková et al. / Computers in Biology and Medicine 43 (2013) 2177–2185

2181

Fig. 6. Expiration (t ¼1 s), left: the velocity (m/s), right: the pressure distribution (Pa).

Table 1 Simulation of breathing in a healthy larynx. Q and R are measured at the bottom boundary Γ bottom (trachea) for inspiration, and at Γ top (larynx–pharynx passage) for expiration in time t¼ 1 s. Breathing process

Max. vel. jvj Flow rate (m/s) Q (l/s)

Inspiration 15.31 Expiration 14.16

1.82 1.66

Pressure Δp (Pa)

Resistance R ðkg=m4 sÞ

Mach no. Ma ½1

106.5 97

0.585 105 0.583 105

4.34 10 2 4.01 10 2

Finally, the Mach number Ma ¼ jvj=cs relates the maximal velocity magnitude jvj to the local speed of sound cs. For the barotropic, adiabatic theﬃ local speed of sound equals pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ pﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ ﬂow cs ¼ γ p=ρ cγργ 1 . In Table 1 we present values for the quantities deﬁned above. This data has been computed for the original geometry (healthy larynx). All data was taken at time t¼1 s for the breathing processes, i.e., when the volumetric ﬂowrate achieves its maximum. Note that due to the different cross-sectional areas on the top and the bottom, the (maximal) inspiration ﬂow rate is larger than the (maximal) expiration ﬂow rate. In our simulations, the variation of density is only of order 10 2 kg=m3 . Thus, in order to compute Reynold's number Re ¼ U n Ln ρn =μn [8], we consider the characteristic density ρn ¼ 1:138 kg=m3 (the density of air) and for the characteristic velocity the maximal inﬂow velocity U n ¼ 4 m=s. This gives Reynold's number about Re 6:3 103 : Here we have computed Ln similarly as in [11] for elliptic pipe, i.e., Ln ¼ 4AEll =U Ell 2:48 cm, where AEll (U Ell ) is the ellipse's surface (circumference, resp.). Due to such high values of Re transient ﬂow up to turbulences can be expected in the ﬂow. We have to note that Reynold's number Re varies signiﬁcantly within the computational domain as well as with respect to time. Particularly, computed locally for the region of vocal folds it is much higher. 3.1. Modiﬁcation of the larynx geometry As described above the aim of this paper is to model the situation where one of the vocal folds is ﬁxed to the speaking

Table 2 Modiﬁcations of the domain. All lengths are in cm. The lengths of the reference post-surgical domain (one vocal fold ﬁxed into the speaking position) are in bold face. T

A

D

T

A

D

0.75 0.95 1.08 1.55

0.65 0.65 0.65 0.65

0.55 0.55 0.55 0.55

0.75 0.95 1.08 1.55

0.65 0.65 0.65 0.65

0.75 0.75 0.75 0.75

0.75 0.95 1.08 1.55

0.55 0.55 0.55 0.55

0.55 0.55 0.55 0.55

0.75 0.95 1.08 1.55

0.55 0.55 0.55 0.55

0.75 0.75 0.75 0.75

0.75 0.95 1.08 1.55

0.47 0.47 0.47 0.47

0.55 0.55 0.55 0.55

0.75 0.95 1.08 1.55

0.47 0.47 0.47 0.47

0.75 0.75 0.75 0.75

position through the injection of fat or another substance. Mathematically, this leads to a change of the computational domain, see Fig. 3. Due to this operation, the size of the gap between the vocal folds is decreased, which affects the air ﬂow to the lungs signiﬁcantly. We perform several simulations of this process and try to identify the most suitable geometry with respect to the quality of breathing. Therefore, we study the air ﬂow in several modiﬁed geometries, evaluate and compare the ﬂow rate Q and the airﬂow resistance R for each modiﬁcation of the domain. Finally, we discuss an optimal design of the larynx shape for the case of ﬁxing one vocal fold into the speaking position. The modiﬁcation of the domain is realized by varying three lengths T, D and A, see Fig. 3 and Section 1. The thickness of the false vocal cords D is always measured at the horizontal coordinate x¼ 2.5. The shortest horizontal distance between the false vocal cords is considered as distance of the false vocal cords T. The length A is measured vertically between the most inner points of the false vocal cords and the healthy vocal folds. Table 2 shows the precise changes of the lengths T; D; A, which have been used in our simulations. We evaluate the data for each domain conﬁguration in the next section.

2182

A. Hundertmark-Zaušková et al. / Computers in Biology and Medicine 43 (2013) 2177–2185

3.2. Data evaluation 3.2.1. Inspiration First we evaluate the inspiration. For this purpose we obtained the airﬂow resistance R and the volumetric ﬂow rate Q for several geometries at time t ¼1 s, see Table 2. We seek a conﬁguration of the larynx geometry that minimizes the resistance and/or maximizes the ﬂow rate. We compare the data computed in modiﬁed geometries with the data obtained in the geometry of a healthy larynx, see Table 1. We started by ﬁxing the thickness of the false vocal cords to D ¼0.55 cm and executed simulations for all chosen distances between the false vocal cords T¼ 0.75, 0.95, 1.08, 1.55 cm and different vertical distances between the false vocal cords and the vocal folds A¼0.47, 0.55, 0.65 cm consecutively (see left column of Table 2). Then we changed D to 0.75 cm and computed the same experiments again (right column of Table 2). The resulting dependencies of the airﬂow resistance R and the ﬂow rate Q on T are shown in Fig. 7. Note that different curves represent different values of A. The reference post-surgical geometry is represented by the curve for A¼0.65 cm at the top of Fig. 7 for T ¼1.08 cm. Additionally, we show the dependence of R and Q on T for both D¼ 0.55 and 0.75 cm for ﬁxed size A¼ 0.55 cm in Fig. 8. The graphical comparisons for A ¼0.47 and A¼ 0.65 cm are quantitatively similar and we omit them here. The ﬁrst consequence we can extract from the graphics in Fig. 8 is that increasing the false vocal cords’ thickness D does not improve the quality of inspiration; the resistance for D¼0.75 cm is higher for all A and the ﬂow rate is smaller. The only exception is the case A¼ 0.47; T¼0.95 cm, cp. Fig. 7 (bottom, right), where the ﬂow rate is slightly higher than for D¼0.55 cm and achieves its local maximum

with respect to T for A¼0.47 cm (Q¼1.827 l/s). This geometry line-up with an increased thickness of the false vocal cords, constricted spacing between the false vocal folds to T¼0.95 cm and a reduced gap between the false vocal cords and the vocal folds (0.47 cm) particularly improves (enlarges) the amount of air in the lungs. The airﬂow resistance, however, achieves in this case 1:95 105 kg=m4 s, which is signiﬁcantly higher than in the case D¼0.55 cm (reference case), cf. Fig. 7 (left). Moreover, we can observe that the highest ﬂow rate is achieved in the case D¼ 0.55, A¼0.55 cm. For a spacing of T ¼1.55 cm between the vocal folds we obtain the maximal ﬂow rate (Q 4 2 l=s), see Fig. 7 (top). However, also the resistance achieves its maximum with respect to T at T ¼1.55 cm. Because of high resistance and since the achieved ﬂow rate exceeds the expected value for a healthy larynx, we exclude the case T ¼1.55 cm from our considerations about the optimal shape. In contrast, comparing the almost constant behavior of the ﬂow rate for T A f0:75; 0:95; 1:08g we can conclude that the optimized larynx geometry for inspiration can be achieved from the reference post-surgical domain with preserving the size A or decreasing it from A ¼0.65 to (e.g.) A¼ 0.55 cm and constricting the distance between the false vocal cords from T ¼1.08 to (e.g.) T¼ 0.75 cm. The thickness D ¼0.55 cm of the false vocal cords remains untouched. These surgical changes of the larynx would minimize the airﬂow resistance and preserve the ﬂow rate about the value in healthy larynx 1:82 l=s (see Table 1).

3.2.2. Expiration The ﬂow rate at Γ bottom during the expiration is dictated by the Dirichlet boundary condition for the velocity. It is not affected by

Fig. 7. Inspiration: dependence of airﬂow resistance (left) and ﬂow rate (right) on the distance of the false vocal cords T for D ¼ 0.55 (top) and D ¼ 0.75 cm (bottom). Case 1: Thickness of false vocal cords D ¼ 0.55 cm. Case 2: Thickness of false vocal cords D¼ 0.75 cm.

A. Hundertmark-Zaušková et al. / Computers in Biology and Medicine 43 (2013) 2177–2185

2183

Fig. 8. Inspiration: dependence of airﬂow resistance (left) and ﬂow rate (right) on the distance of false vocal cords T for A ¼ 0.55 cm.

Fig. 9. Expiration: dependence of the airﬂow resistance on the distance of false vocal cords T for D ¼ 0.55 (left) and for D¼ 0.75 cm (right). (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this paper.)

modiﬁcation of the lengths T, D, A, since it is measured before the air passes the modiﬁed part of the larynx, hence, it is not relevant for our study. However, we have observed that the ﬂow rates measured at the top boundary Γ top are also only nonsigniﬁcanty affected by geometry modiﬁcation: the values vary in the range ½1:655; 1:663 l=s and achieve, therefore, values of expiration ﬂow rate from healthy larynx, cp. Table 1. Due to the above reasons the maximization of the ﬂow rate fails off and we omit its graphical presentation as a function of T in what follows. Consequently, the resistance during the expiration process is governed mainly by the pressure difference. In Fig. 9 we can observe clearly that the lowest resistance for the gap size A¼ 0.65 cm is obtained when the distance of the vocal folds is T¼ 1.08 cm, i.e., in the reference post-surgical geometry (the green line with circles for T¼1.08). We also observe a local resistance minimum for the distance T¼ 0.95 cm when the gap size A¼ 0.55 cm and the thickness D¼0.75 cm, though the difference in resistance for other gap sizes and/or thickness D¼0.55 cm is of low magnitude. Furthermore, this resistance is only slightly higher than the value in the reference post-surgical geometry, see Fig. 9. Comparing the results in Fig. 9 we conclude that changing the thickness D does not play a signiﬁcant role in improving the process of expiration. We close this section with an overview of our study's results on the volumetric ﬂow rate and the airﬂow resistance, see Fig. 10. Each marker represents one geometry. We denoted the different lengths T by four different colors. The shapes represent different values for A, i.e., square for A¼0.47 cm; triangle for A¼0.55 cm and circle for A¼ 0.65 cm. The blank markers represent D¼ 0.55 cm, the ﬁlled ones represent D¼0.75 cm.

4. Discussion As we can deduce from Fig. 10 we do not ﬁnd a unique optimal geometry for reducing the resistance and increasing the ﬂow rate at once. Instead, we obtain several geometries that yield either maximal ﬂow rate or minimal resistance, though, the ﬂow rate remains almost constant throughout the study. Reducing T to 0.75 and A to 0.55 cm leads to a decrease of the airﬂow resistance for inspiration from R ¼ 1:940 105 in the reference post-surgical geometry to R ¼ 1:807 105 kg=m4 s while the ﬂow rate stays about the same as in the reference post-surgical domain or in the healthy larynx. We obtain similar results by keeping A at 0.65 cm and changing T to 0.75 cm as above, see Fig. 7. Generally, we have observed that increasing the thickness of the false vocal folds to D¼ 0.75 cm increases the airﬂow resistance for inspiration. When we, additionally, decrease T to 0.95 cm and reduce the gap size to A ¼0.47 cm we obtain an increase of the volumetric ﬂow rate during inspiration of about 10 ml/s in comparison to the reference case. This even exceeds the value for the healthy larynx. However, because of the high resistance (even higher than in the reference post-surgical case) the signiﬁcance of the ﬂow rate increase by this modiﬁcation of the geometry is disputable. Taking into account that the resistance R increases almost linearly with respect to the vestibular folds’ distance T while the ﬂow rate Q remains almost constant for T A f0:75; 0:95; 1:08g (see Fig. 7) we arrive at an optimal larynx shape for the inspiration process: for T¼ 0.75, A¼ 0.55 (0.65), D¼0.55 cm the ﬂow rate Q¼1.815 (1.806) l/s is close to the case of the healthy larynx and the resistance R ¼ 1:81ð1:79Þ 105 kg=m4 s reaches the study's minimum.

2184

A. Hundertmark-Zaušková et al. / Computers in Biology and Medicine 43 (2013) 2177–2185

Fig. 10. Overview on ﬂow rates and resistances for all tested larynx geometries (measured at Γ bottom for inspiration and at Γ top for expiration). Top: inspiration versus expiration, bottom: details on inspiration (left) and expiration (right). (For interpretation of the references to color in this ﬁgure caption, the reader is referred to the web version of this paper.

Focusing on an optimal larynx shape for the expiration process the reference post-surgical geometry (T ¼1.08, A ¼0.65, D ¼0.55 cm) provides minimal resistance R ¼ 1:291 105 kg=m4 s whereas the ﬂow rate Q achieves the same value as for the healthy larynx, Q¼1.656 l/s. Another local resistance minimum is obtained at R ¼ 1:338 105 kg=m4 s by reducing T to 0.95 cm and A to 0.47 cm. Increasing the thickness D to 0.75 cm does not have a signiﬁcant inﬂuence on the expiration airﬂow resistance. However, we cannot conclude that one of these geometries is optimal for the inspiration process as well. The inspiration resistance in these cases is signiﬁcantly higher compared to other geometries. Now we want to derive an improved geometry for the whole breathing process from the obtained data balancing both inspiration and expiration. Since the expiration resistances for T ¼0.95 and T¼ 1.08 cm do not differ signiﬁcantly (cp. Fig. 9) we chose T ¼0.95 cm as compromise. For the inspiration process this still yields the second lowest value for the resistance, see Figs. 7 and 8,

without strongly affecting the ﬂow rate. Reducing the gap to A¼0.55 cm and increasing the thickness to D¼ 0.75 cm would provide almost optimal expiration resistance, however, the drawback for the inspiration is relatively big. Therefore, we arrive at the optimal geometry T ¼0.95, A ¼0.65, D¼0.55 cm. We summarize the initial (healthy larynx), the reference postsurgical and the ﬁnal optimized geometry for inspiration and expiration.

5. Conclusion The goal of this paper is to present results on mathematical modeling and numerical simulation of glottal ﬂow. We simulate and qualify the airﬂow in a human larynx. The simpliﬁed airﬂow model consisting of the weakly compressible Navier–Stokes equations in the two dimensional domain obtained after vertical larynx

A. Hundertmark-Zaušková et al. / Computers in Biology and Medicine 43 (2013) 2177–2185

Acknowledgments

Table 3 Optimized geometries for breathing (inspiration and expiration). Geometry

Inspiration

Expiration

Q ðl=sÞ

R ðkg=m4 sÞ

Q ðl=sÞ

R ðkg=m4 sÞ

0.585 105

1.655

0.583 105

1.940 105

1.656

1.291 105

1.871 105 1.917 105

1.660 1.657

1.371 105 1.334 105

1.656

1.291 105

Healthy larynx 1.820 Post-surgical reference T ¼ 1:08; A ¼ 0:65; D ¼ 0:55 1.817 Optimized for both inh. and exh. T ¼ 0:95; A ¼ 0:65; D ¼ 0:55 1.814 T ¼ 0:95; A ¼ 0:55; D ¼ 0:75 1.811 Optimal for inspiration T ¼ 0:75; A ¼ 0:55; D ¼ 0:55 1.815 T ¼ 0:75; A ¼ 0:65; D ¼ 0:55 1.806 Optimal for expiration T ¼ 1:08; A ¼ 0:65; D ¼ 0:55

1.807 105 1.791 105

cut has been used. The ﬁnite element method implemented in COMSOL has been applied for numerical simulations. The quality of breathing in the human larynx has been quantiﬁed using the volumetric ﬂow rate Q and the airﬂow resistance R. These values have been measured at ﬁrst for the case of a healthy larynx. Then we designed a reference post-surgical geometry with an enlarged left vocal fold. For this geometry, particularly the airﬂow resistance increased, see Table 3. Hence, we modiﬁed the post-surgical domain by varying the (horizontal) distance T between the false vocal cords, the thickness D of the false vocal cords and the size of the (vertical) gap A between the vocal folds and the false vocal folds to maximize the ﬂow rate and to minimize the airﬂow resistance for inspiration and expiration. Our numerical simulations did not deliver a unique optimal geometry for both inspiration and expiration at once. The optimal geometry for inspiration does not improve the quality of expiration compared to our reference post-surgical geometry and vice versa, whereas the airﬂow in the reference geometry during the expiration shows the optimal behavior. Nevertheless, we can improve the quality of inspiration, i.e., decrease the airﬂow resistance considerably by reducing the distance of the false vocal folds from 1.08 to 0.95 cm. This geometry modiﬁcation only slightly increases the expiration airﬂow resistance compared to the optimal geometry. Therefore, the geometry setup T¼ 0.95, A¼0.65, D¼0.55 cm balances best the resistance and airﬂow for inspiration and expiration. In order to approximate the aerodynamics in the laryngeal ﬂow more precisely, we point out possible improvements of its mathematical and numerical modeling, which will be a goal of our future study. Firstly, one can consider more realistic MRI-scan-based three dimensional geometries. Physiologically dictated pressure perturbation boundary condition may be considered. Moreover, boundary conditions originating from the coupling of 3D with 1D/0D model would take into account the remaining part of the respiratory system. Furthermore, techniques of mathematical optimization can be used in order to recommend a new, surgically changed larynx shape. Eventually, further aspects and speciﬁc patient requirements should be reﬂected. Conﬂict of interest statement None declared.

2185

The authors thank Mária Lukáčová-Medvid’ová from the University of Mainz for mathematical discussions and enthusiastic support and also Irina Brucker and Alexandra Ickert for their help by working on our manuscript. The present research has been partially supported under the DFG project ZA 613/1-1.

References [1] F. Alipour, R.C. Scherer, Pulsatile airﬂow during phonation: an excised larynx model, J. Acoust. Soc. Am. 97 (2) (1995) 1241–1248. [2] F. Alipour, R.C. Scherer, Flow separation an a computational oscillating vocal fold model, J. Acoust. Soc. Am. 116 (3) (2004) 1710–1719. [3] F. Alipour, R.C. Scherer, Ventricular pressures in phonating excised larynges, J. Acoust. Soc. Am. 132 (2) (2012) 1017–1726. [4] A.G. Apostoli. Dynamic Modelling of the Human Larynx in Phonation, Ph.D. Thesis, University of Edingburgh, 2012. [5] L. Bailly, X. Pelorson, N. Henrich, N. Ruty, Inﬂuence of a constriction in the near ﬁeld of the vocal folds: physical modeling and experimental validation, J. Acoust. Soc. Am. 124 (5) (2008) 3296–3308. [6] N.E. Chisari, G. Artana, D. Sciamarella, Experimental and numerical study of patterns in laryngeal ﬂow, J. Phys. Conf. Ser. 166 (2009) 012013. [7] COMSOL Multhiphysics, 〈http://www.comsol.com/products/multiphysics/〉. [8] M. Feistauer, J. Felcman, I. Straškraba, Mathematical and Computational Methods for Compressible Flow, Oxford Science Publications, 2003. [9] H. Gray, Anatomy of the Human Body, 20th ed., 1918. [10] John Hopkins School of Medicine's Interactive Respiratory Physiology. Airway Resistance 〈http://oac.med.jhmi.edu/res_phys/Encyclopedia/AirwayResis tance/AirwayResistance.HTML〉. [11] S.Y. Huang, F. Mayinger, Wärmeübergang bei freier Konvektion um elliptische Rohre, in: Wärme- und Stoffübertragung, vol. 18, Springer Verlag, 1984, pp. 175–183. [12] L. Int Panis, B. de Geus, G. Vandenbulcke, H. Willems, B. Degraeuwe, N. Bleux, Exposure to UFP and PM in trafﬁc: a comparison of cyclists and car drivers, Atmos. Environ. 44 (2010) 2263–2270. [13] K. Kadoya, N. Matsunaga, A. Nagashima, Viscosity and thermal conductivity of dry air in the gaseous phase, J. Phys. Chem. Ref. Data 14 (4) (1985) 947–970. [14] S. Khosla, S. Muruguppan, E. Gutmark, R. Scherer, Vortical ﬂow ﬁeld during phonation in an excised canine larynx model, Ann. Otol. Rhinol. Laryngol. 116 (3) (2007) 217–228. [15] R. Langnickel, Experimentelle Untersuchungen zur Frage des Strömungswiderstandes am Larynx-Trachealpräparat bei normaler und veränderter Stimmritze, Eur. Arch. Oto-Rhino-Laryngol. 212 (1) (1976) 43–56. [16] S. Li, R.C. Scherer, M. Wan, S. Wang, The effect of entrance radii on intraglottal pressure distributions in the divergent glottis, J. Acoust. Soc. Am. 131 (2) (2012) 1371–1377. [17] H. Luo, X. Zheng, R. Mittal, H. Dong, S. Bielamowicz, R. Walsh, J. Hahn, An immersed-boundary method for ﬂow–structure interaction in biological systems with application to phonation, J. Comput. Phys. 227 (22) (2008) 9303–9332. [18] H. Nomura, T. Funada, Numerical simulation of unsteady ﬂow through the rigid glottis, Acoust. Sci. Technol. 27 (3) (2006). [19] X. Pelorson, A. Hirschberg, A.P. Wijnands, H.Bailliet. Description of the ﬂow through in-vitro models of the glottis during phonation, Acta Acoust. 3 (1995) 191–202. [20] J. Prokopová, M. Feistauer, J. Horáček, V. Kučera, Numerical simulation of airﬂow through the model of oscillating vocal folds, Appl. Comput. Mech. 4 (2010) 215–224. [21] R.C. Scherer, D. Shinwari, K.J. De Witt, C. Zhang, B.R. Kucinschi, A.A. Afjeh, Intraglottal pressure proﬁles for a symmetric and oblique glottis with a divergence angle of 10 degrees, J. Acoust. Soc. Am. 109 (4) (2001) 1616–1630. [22] R.C. Scherer, K.J. De Witt, B.R. Kucinschi, The effect of exit radii on intraglottal pressure distributions in the convergent glottis, J. Acoust. Soc. Am. 110 (5 Pt 1) (2001) 2267–2269. [23] B.H. Story, I.R. Titze, Voice simulation with a body-cover model of the vocal folds, J. Acoust. Soc. Am. 97 (2) (1995) 1249–1260. [24] I.R. Titze, S.S. Schmidt, M.R. Titze, Phonation threshold pressure in a physical model of the vocal fold mucisa, J. Acoust. Soc. Am. 97 (5) (1995) 3080–3084.