Computers in Biology and Medicine 43 (2013) 2177–2185
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Numerical simulation of glottal flow 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 difficulties with breathing but rarely with voicing. However, clinical experience shows that the shape of the larynx (voice box) seems to have a significant influence on the degree of airflow and breathing pattern. In order to find 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 flow for weakly compressible Navier–Stokes equations has been performed. We compare airflow resistance and volumetric flow 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 flow Numerical simulation Weakly compressible flow Geometry optimization Respiration Airflow resistance Volumetric flow rate
1. Respiration process and speaking In the respiration process (inspiration and expiration) and during speaking, air flows 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 configuration 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 airflow). For voice production (i.e., phonation) the vocal folds are brought together actively by muscular actions. This configuration of adducted vocal folds enables air pressure to build up under the closed (adducted) vocal folds until they are pushed apart due to their floppy, 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 flowing 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 flow regimes, flow 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 confirmed. 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 significantly affects the pressure distributions on the glottis and needs to be well specified 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 fluid– structure interaction problem is solved. The flow 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 fluid. Langnickel [15] belongs to the first studies about the influence of the changed vocal folds’ shape on the airflow and its resistance. In Fig. 3 a possible configuration 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 modified 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 difficulties (shortness of breath) causing the feeling of lack of air in the lungs. If vocal folds do not actively move (separation or adduction), the ‘fixed’ gap between both vocal folds determines how breathing and voicing is. The more they are separated, the better the condition for respiratory airflow. 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 flow during the respiration process in vitro and to design a suitable geometry with respect to the air volume and flow resistance after passing the modified vocal folds. Here, we assume the left vocal fold to be fixed 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 simplified 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 configurations. 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 fixed 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 simplified 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 flow can be modeled using the weakly compressible Navier– Stokes equations. The additional assumption of adiabatic and barotropic flow 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 fluid ρ 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 coefficients 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 simplified 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 fitting 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 outflow boundary which is now Γ top we fix 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 finite element method allows to use several types of elements (Lagrange polynomial of the first, second and third orders), for instance the P1–P1, P2–P1 or P3–P2 elements. The presented results are obtained using first order Lagrange polynomial for pressure as well as velocity (P1–P1). For time discretization the first 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 simplified two dimensional cut, see Fig. 4 (left). Figs. 5 and 6 show the streamlines, the velocity field 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 modified larynx geometries. In what follows we intend to find the optimal geometry with respect to these quantities. In order to measure the volume of air flowing into the lungs, R we use the flow 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 outflow boundary Γ bottom equals the transverse elliptic diameter a~ ¼ 2a ¼ 2:5 cm, thus, we obtain an approximation of the volumetric flow rate Q: Z v nΓ bottom : Q b~ Q~ ¼ b~ Γ bottom
The flow rate is computed in liters per second. During the expiration it is useful to evaluate the flow rate above the vocal cords, i.e., the developed flow on the top boundary Γ top . Note that the Dirichlet boundary condition for velocity at the bottom strongly affects the flow 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 define the airflow 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 flow rate, cf. [10,15]. For optimization we will focus on the flow rate Q and the airflow 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 theffi local speed of sound equals pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi flow cs ¼ γ p=ρ cγργ 1 . In Table 1 we present values for the quantities defined 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 flowrate achieves its maximum. Note that due to the different cross-sectional areas on the top and the bottom, the (maximal) inspiration flow rate is larger than the (maximal) expiration flow 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 inflow 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 flow up to turbulences can be expected in the flow. We have to note that Reynold's number Re varies significantly 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. Modification of the larynx geometry As described above the aim of this paper is to model the situation where one of the vocal folds is fixed to the speaking
Table 2 Modifications of the domain. All lengths are in cm. The lengths of the reference post-surgical domain (one vocal fold fixed 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 flow to the lungs significantly. 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 flow in several modified geometries, evaluate and compare the flow rate Q and the airflow resistance R for each modification of the domain. Finally, we discuss an optimal design of the larynx shape for the case of fixing one vocal fold into the speaking position. The modification 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 configuration 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 airflow resistance R and the volumetric flow rate Q for several geometries at time t ¼1 s, see Table 2. We seek a configuration of the larynx geometry that minimizes the resistance and/or maximizes the flow rate. We compare the data computed in modified geometries with the data obtained in the geometry of a healthy larynx, see Table 1. We started by fixing 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 airflow resistance R and the flow 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 fixed 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 first 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 flow rate is smaller. The only exception is the case A¼ 0.47; T¼0.95 cm, cp. Fig. 7 (bottom, right), where the flow 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 airflow resistance, however, achieves in this case 1:95 105 kg=m4 s, which is significantly higher than in the case D¼0.55 cm (reference case), cf. Fig. 7 (left). Moreover, we can observe that the highest flow 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 flow 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 flow 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 flow 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 airflow resistance and preserve the flow rate about the value in healthy larynx 1:82 l=s (see Table 1).
3.2.2. Expiration The flow 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 airflow resistance (left) and flow 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 airflow resistance (left) and flow rate (right) on the distance of false vocal cords T for A ¼ 0.55 cm.
Fig. 9. Expiration: dependence of the airflow 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 figure caption, the reader is referred to the web version of this paper.)
modification of the lengths T, D, A, since it is measured before the air passes the modified part of the larynx, hence, it is not relevant for our study. However, we have observed that the flow rates measured at the top boundary Γ top are also only nonsignificanty affected by geometry modification: the values vary in the range ½1:655; 1:663 l=s and achieve, therefore, values of expiration flow rate from healthy larynx, cp. Table 1. Due to the above reasons the maximization of the flow 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 significant role in improving the process of expiration. We close this section with an overview of our study's results on the volumetric flow rate and the airflow 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 filled ones represent D¼0.75 cm.
4. Discussion As we can deduce from Fig. 10 we do not find a unique optimal geometry for reducing the resistance and increasing the flow rate at once. Instead, we obtain several geometries that yield either maximal flow rate or minimal resistance, though, the flow rate remains almost constant throughout the study. Reducing T to 0.75 and A to 0.55 cm leads to a decrease of the airflow resistance for inspiration from R ¼ 1:940 105 in the reference post-surgical geometry to R ¼ 1:807 105 kg=m4 s while the flow 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 airflow 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 flow 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 significance of the flow rate increase by this modification 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 flow 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 flow 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 flow 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 figure 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 flow 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 significant influence on the expiration airflow 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 significantly 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 significantly (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 flow 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 final 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 flow. We simulate and qualify the airflow in a human larynx. The simplified airflow 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 finite element method implemented in COMSOL has been applied for numerical simulations. The quality of breathing in the human larynx has been quantified using the volumetric flow rate Q and the airflow resistance R. These values have been measured at first 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 airflow resistance increased, see Table 3. Hence, we modified 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 flow rate and to minimize the airflow 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 airflow in the reference geometry during the expiration shows the optimal behavior. Nevertheless, we can improve the quality of inspiration, i.e., decrease the airflow resistance considerably by reducing the distance of the false vocal folds from 1.08 to 0.95 cm. This geometry modification only slightly increases the expiration airflow 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 airflow for inspiration and expiration. In order to approximate the aerodynamics in the laryngeal flow 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 specific patient requirements should be reflected. Conflict 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 airflow 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, Influence of a constriction in the near field 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 flow, 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 traffic: 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 flow field 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 flow–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 flow through the rigid glottis, Acoust. Sci. Technol. 27 (3) (2006). [19] X. Pelorson, A. Hirschberg, A.P. Wijnands, H.Bailliet. Description of the flow 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 airflow 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 profiles 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.
Contents lists available at ScienceDirect
Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm
Numerical simulation of glottal flow 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 difficulties with breathing but rarely with voicing. However, clinical experience shows that the shape of the larynx (voice box) seems to have a significant influence on the degree of airflow and breathing pattern. In order to find 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 flow for weakly compressible Navier–Stokes equations has been performed. We compare airflow resistance and volumetric flow 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 flow Numerical simulation Weakly compressible flow Geometry optimization Respiration Airflow resistance Volumetric flow rate
1. Respiration process and speaking In the respiration process (inspiration and expiration) and during speaking, air flows 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 configuration 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 airflow). For voice production (i.e., phonation) the vocal folds are brought together actively by muscular actions. This configuration of adducted vocal folds enables air pressure to build up under the closed (adducted) vocal folds until they are pushed apart due to their floppy, 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 flowing 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 flow regimes, flow 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 confirmed. 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 significantly affects the pressure distributions on the glottis and needs to be well specified 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 fluid– structure interaction problem is solved. The flow 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 fluid. Langnickel [15] belongs to the first studies about the influence of the changed vocal folds’ shape on the airflow and its resistance. In Fig. 3 a possible configuration 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 modified 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 difficulties (shortness of breath) causing the feeling of lack of air in the lungs. If vocal folds do not actively move (separation or adduction), the ‘fixed’ gap between both vocal folds determines how breathing and voicing is. The more they are separated, the better the condition for respiratory airflow. 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 flow during the respiration process in vitro and to design a suitable geometry with respect to the air volume and flow resistance after passing the modified vocal folds. Here, we assume the left vocal fold to be fixed 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 simplified 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 configurations. 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 fixed 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 simplified 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 flow can be modeled using the weakly compressible Navier– Stokes equations. The additional assumption of adiabatic and barotropic flow 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 fluid ρ 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 coefficients 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 simplified 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 fitting 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 outflow boundary which is now Γ top we fix 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 finite element method allows to use several types of elements (Lagrange polynomial of the first, second and third orders), for instance the P1–P1, P2–P1 or P3–P2 elements. The presented results are obtained using first order Lagrange polynomial for pressure as well as velocity (P1–P1). For time discretization the first 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 simplified two dimensional cut, see Fig. 4 (left). Figs. 5 and 6 show the streamlines, the velocity field 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 modified larynx geometries. In what follows we intend to find the optimal geometry with respect to these quantities. In order to measure the volume of air flowing into the lungs, R we use the flow 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 outflow boundary Γ bottom equals the transverse elliptic diameter a~ ¼ 2a ¼ 2:5 cm, thus, we obtain an approximation of the volumetric flow rate Q: Z v nΓ bottom : Q b~ Q~ ¼ b~ Γ bottom
The flow rate is computed in liters per second. During the expiration it is useful to evaluate the flow rate above the vocal cords, i.e., the developed flow on the top boundary Γ top . Note that the Dirichlet boundary condition for velocity at the bottom strongly affects the flow 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 define the airflow 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 flow rate, cf. [10,15]. For optimization we will focus on the flow rate Q and the airflow 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 theffi local speed of sound equals pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffi flow cs ¼ γ p=ρ cγργ 1 . In Table 1 we present values for the quantities defined 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 flowrate achieves its maximum. Note that due to the different cross-sectional areas on the top and the bottom, the (maximal) inspiration flow rate is larger than the (maximal) expiration flow 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 inflow 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 flow up to turbulences can be expected in the flow. We have to note that Reynold's number Re varies significantly 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. Modification of the larynx geometry As described above the aim of this paper is to model the situation where one of the vocal folds is fixed to the speaking
Table 2 Modifications of the domain. All lengths are in cm. The lengths of the reference post-surgical domain (one vocal fold fixed 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 flow to the lungs significantly. 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 flow in several modified geometries, evaluate and compare the flow rate Q and the airflow resistance R for each modification of the domain. Finally, we discuss an optimal design of the larynx shape for the case of fixing one vocal fold into the speaking position. The modification 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 configuration 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 airflow resistance R and the volumetric flow rate Q for several geometries at time t ¼1 s, see Table 2. We seek a configuration of the larynx geometry that minimizes the resistance and/or maximizes the flow rate. We compare the data computed in modified geometries with the data obtained in the geometry of a healthy larynx, see Table 1. We started by fixing 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 airflow resistance R and the flow 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 fixed 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 first 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 flow rate is smaller. The only exception is the case A¼ 0.47; T¼0.95 cm, cp. Fig. 7 (bottom, right), where the flow 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 airflow resistance, however, achieves in this case 1:95 105 kg=m4 s, which is significantly higher than in the case D¼0.55 cm (reference case), cf. Fig. 7 (left). Moreover, we can observe that the highest flow 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 flow 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 flow 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 flow 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 airflow resistance and preserve the flow rate about the value in healthy larynx 1:82 l=s (see Table 1).
3.2.2. Expiration The flow 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 airflow resistance (left) and flow 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 airflow resistance (left) and flow rate (right) on the distance of false vocal cords T for A ¼ 0.55 cm.
Fig. 9. Expiration: dependence of the airflow 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 figure caption, the reader is referred to the web version of this paper.)
modification of the lengths T, D, A, since it is measured before the air passes the modified part of the larynx, hence, it is not relevant for our study. However, we have observed that the flow rates measured at the top boundary Γ top are also only nonsignificanty affected by geometry modification: the values vary in the range ½1:655; 1:663 l=s and achieve, therefore, values of expiration flow rate from healthy larynx, cp. Table 1. Due to the above reasons the maximization of the flow 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 significant role in improving the process of expiration. We close this section with an overview of our study's results on the volumetric flow rate and the airflow 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 filled ones represent D¼0.75 cm.
4. Discussion As we can deduce from Fig. 10 we do not find a unique optimal geometry for reducing the resistance and increasing the flow rate at once. Instead, we obtain several geometries that yield either maximal flow rate or minimal resistance, though, the flow rate remains almost constant throughout the study. Reducing T to 0.75 and A to 0.55 cm leads to a decrease of the airflow resistance for inspiration from R ¼ 1:940 105 in the reference post-surgical geometry to R ¼ 1:807 105 kg=m4 s while the flow 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 airflow 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 flow 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 significance of the flow rate increase by this modification 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 flow 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 flow 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 flow 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 figure 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 flow 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 significant influence on the expiration airflow 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 significantly 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 significantly (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 flow 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 final 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 flow. We simulate and qualify the airflow in a human larynx. The simplified airflow 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 finite element method implemented in COMSOL has been applied for numerical simulations. The quality of breathing in the human larynx has been quantified using the volumetric flow rate Q and the airflow resistance R. These values have been measured at first 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 airflow resistance increased, see Table 3. Hence, we modified 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 flow rate and to minimize the airflow 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 airflow in the reference geometry during the expiration shows the optimal behavior. Nevertheless, we can improve the quality of inspiration, i.e., decrease the airflow resistance considerably by reducing the distance of the false vocal folds from 1.08 to 0.95 cm. This geometry modification only slightly increases the expiration airflow 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 airflow for inspiration and expiration. In order to approximate the aerodynamics in the laryngeal flow 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 specific patient requirements should be reflected. Conflict 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 airflow 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, Influence of a constriction in the near field 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 flow, 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 traffic: 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 flow field 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 flow–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 flow through the rigid glottis, Acoust. Sci. Technol. 27 (3) (2006). [19] X. Pelorson, A. Hirschberg, A.P. Wijnands, H.Bailliet. Description of the flow 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 airflow 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 profiles 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.
