Effects of tracheal stenosis on flow dynamics in upper human airways

(1)

c

TU Delft, The Netherlands, 2006

EFFECTS OF TRACHEAL STENOSIS ON FLOW

DYNAMICS IN UPPER HUMAN AIRWAYS

Santhosh T Jayaraju∗, Mark Brouns†, Chris Lacor†, Johan De Mey†† and Sylvia Verbanck††

∗_{†Vrije Universiteit Brussel, Department of Mechanical Engineering}

Pleinlaan 2, 1050 Brussel, Belgium

e-mail: santhosh@stro.vub.ac.be, mark.brouns@vub.ac.be, chris.lacor@vub.ac.be web page: http://stro.vub.ac.be

††_{Academic hospital, Vrije Universiteit Brussel, Laarbeeklaan 101, 1090 Brussel, Belgium}

e-mail: johan.demey@az.vub.ac.be, sylvia.verbanck@az.vub.ac.be

Key words: Human oral airway, Tracheal stenosis, Larynx, RANS

Abstract. In order to investigate the tracheal pressures for different percentages of stenotic constriction, flow simulations are performed on a realistic geometry based on CT-scans of a patient for sedentary breathing (15 l/min) and normal breathing (30 l/min) conditions using a Reynolds Averaged Navier Stokes approach for unstructured hexahedral meshes. Turbulence modeling is based on low Reynolds number Yang-Shih k − ε model. The realistic airway model in itself is highly irregular and the presence of stenosis adds to its geometrical complexity, resulting in very complex flow patterns with flow separations and skewed velocity profiles. Detailed analysis of these flows is presented for 0%, 49%, 75%, 84% and 91% stenotic constrictions. The pressure drop shows modest increases with the degree of narrowing up to 75% constriction, beyond which it steeply rises. Such a pattern agrees with the appearance of breathing symptoms only when patients already show a very marked stenosis. In order to temporarily ease work of breathing, a decrease of the pressure drop is sometimes induced by having the patient inhale heliox. An additional flow analysis is performed with heliox as the working fluid.

1 INTRODUCTION

(2)

Santhosh T Jayaraju, Mark Brouns, Chris Lacor and Sylvia Verbanck

patients post-intubation. The causes of adult laryngeal and upper tracheal stenosis may be trauma, chronic inflammatory diseases or post-operative change1_{. In the present work,} we first concentrate on providing detailed description of the flow dynamics in the realistic geometry without stenosis, followed by analysis in the presence of four different stenotic constrictions.

Figure 1: Schematic representation of airway and the realistic airway geometry

(3)

and noticed that the difference in pressure drop is more dominated by the cross-sectional area than by the shape of the glottis.

A few experimental studies have also been performed recently. Corcoran et al.9 stud-ied the laryngeal jet from the cadaver of an adult female and noted the importance of downstream conditions on the position of laryngeal jet. PIV measurements on a sim-plified model of the throat were performed by Heenan et al.2 _{and compared with CFD} predictions to conclude that the increased viscous effects at lower Reynolds numbers are not adequately captured owing to the deficiencies in CFD modeling. However, very less attention is paid to the study of stenosis. Wassermann et al.10 _{used a newly developed} bronchoscopic technique for the assessment of intratracheal pressures as a way to quantify tracheal resistance for use in diagnosis of patients with tracheal stenosis. Fasano et al.11 studied the frequency dependence of specific airway resistance in patients with laryngeal obstructions. To the author’ss knowledge, the present work is the first attempt using CFD to provide comprehensive characterization of airflow in realistic geometries with stenosis.

2 MODEL PREPARATION

Figure 2: Sample of a CT-scan

(4)

Figure 3: View of unstructured hexahedral mesh

The highly irregular and complex nature of the model makes the creation of structured grid tedious and work-intensive. For such complex geometries, unstructured meshes are more suited. Using Numeca’s unstructured grid generator - Hexpress (2.2β), an all-hexahedral unstructured mesh was generated, containing approximately 750,000 cells. Fig. 3 shows the surface mesh and a zoom in the region of mouth-throat. Box-adaptation was used to locally refine the mesh near the tracheal stenosis.

3 NUMERICAL METHODS

The computational domain is imported into Numeca (Fine/Hexa 2.1β) compressible RANS solver. The basic Reynolds-Averaged-Navier-Stokes equation in a cartesian frame of reference integrated over a control volume V is expressed as,

Z Z V Z _∂U ∂tdV + Z Z S ~ F · ~dS = Z Z S ~ Q · ~dS (1)

(5)

The flow is modeled using a Low-Reynolds number Yang-Shih k − ε turbulence model. The accuracy of the model is tested by performing a simulation on idealized mouth-throat geometry of Heenan et al.2 _{for which experimental results are available. For a flow rate} of 30 l/min, the measured and predicted pressure drops agreed well (Measured: 30 Pa, Simulation: 29.2 Pa) which indicates that the selected turbulence model is appropriate. More details on the turbulence model can be found in the Numeca’s user manual17_{. The} value of turbulent kinetic energy at the inlet is prescribed assuming 5% turbulence inten-sity. Kleinstreuer et al.18 _{studied the effects of turbulence intensities (in the range of 1} -10%) on the flow in an idealized mouth-throat geometry constructed based on the data set of Cheng et al.19_{. It was seen that the variation of turbulence intensity had negligible} effect of pressure drop along the domain for a 30 l/min flow rate.

4 RESULTS

(6)

(a) (b)

(c) (d) (e)

(f) (g) (h)

(7)

(a) (b)

(c) (d) (e)

(f) (g) (h)

(8)

(a) (b)

(c) (d) (e)

(f) (g) (h)

Figure 6: Normalized kinetic energy profiles (k/u2

in) in the oral airway model. (a) At 15 l/min. (b) At

(9)

a transition to turbulence in the downstream regions of the shear layer which modifies the location of reattachment. Also, the secondary motions at sections E1-E2 and F1-F2 reduce when compared to that of 15 l/min.

Fig. 6.a shows the normalized turbulent kinetic energy along the domain for 15 l/min. The maximum value is at the inlet (prescribed assuming 5% turbulence intensity) which is damped out before the flow reaches half way through the mouth and laminar flow prevails through the rest of the domain. Similar behavior is reported by Kleinstreuer et al.18_{. At a flow rate of 30 l/min, the turbulent intensity magnifies soon after the glottis} region, reaching a maximum value in between 1-3 diameters downstream of glottis. The maximum turbulence intensity is concentrated on the posterior side at one diameter from the glottis (D1-D2), moves towards the centre at three diameters (E1-E2) and finally towards the anterior side at six diameters downstream (F1-F2). These results are in contrast with the observations of Kleinstreuer et al.18 _{and Corcoran et al.}9 _{who report} maximum fluctuations in the shear layers close to the walls. This might be due to over-simplification of glottis and tracheal region which changes the evolution of the laryngeal jet when compared to the realistic geometry.

(a)

Figure 7: Area-averaged pressure drop (p − pin) taken at every 5mm along the oral airway model

(10)

comparatively much less and remains constant through the trachea as the laminar flow prevails.

(a) 49% Stenosis (b) 75% Stenosis (c) 84% Stenosis (d) 91% Stenosis

(e) (f) (g) (h)

(i) (j) (k) (l)

Figure 8: Velocity profiles for a flow rate of 15 l/min. (a)-(d) 2D streaklines. (e)-(h) Axial velocity contour (magnitudes in cm/s) and secondary velocity vector lines at A1-A2. (i)-(l) Axial velocity contour and secondary velocity vector lines at B1-B2

(11)

the flow separation is much bigger (Fig. 8.b) and also the secondary vortex is more profound even three diameters downstream (Fig. 8.j). The turbulence intensity for both 49% and 75% constriction remains very small (Fig. 10.a) showing that there is no onset of turbulence up to 75% stenosis. However, the flow turns turbulent for 84% and 91% with a fully developed recirculation region on the posterior side (Fig. 8.c,d) with re-attachment lengths of about three diameters downstream. From Fig. 8.(e)-(h) we see that the magnitude of velocity on the anterior side increases with increase in constriction percentage. With 84% and 91% constriction, we also see a pocket of high velocity region on the posterior side (Fig. 8.g,h) which is due to stronger flow recirculation.

(a) 49% Stenosis (b) 75% Stenosis (c) 84% Stenosis (d) 91% Stenosis

(e) (f) (g) (h)

(i) (j) (k) (l)

Figure 9: Velocity profiles for a flow rate of 30 l/min. (a)-(d) 2D streaklines. (e)-(h) Axial velocity contour (magnitudes in cm/s) and secondary velocity vector lines at A1-A2. (i)-(l) Axial velocity contour and secondary velocity vector lines at B1-B2

(12)

the flow separation region on the posterior wall for 49% and 75% constriction (Fig. 9.a,b) are quite different when compared to those of 15 l/min (Fig. 8.a,b). This difference is due to the fact that the flow is already turbulent for 30 l/min whereas it still remains laminar for 15 l/min. For 84% and 91% constrictions, the flow is turbulent for both flow rates and hence we observe similar structures. However, the recirculation region (Fig. 9.c,d) is slightly bigger with delayed flow re-attachment when compared to those of 15 l/min (Fig. 8.c,d) which is understandably due to higher mass flow rate. The same is the reason for more profound secondary velocities at three diameters downstream (Fig. 9.k,l) when compared to those of 15 l/min (Fig. 8.k,l).

(a) (b)

Figure 10: Area-averaged normalized kinetic energy (k/u2

in) at every 5mm after the stenotic constriction.

(a) 15 l/min (b) 30 l/min

In Fig. 10.b we see that the turbulence intensity for 49% constriction remains slightly higher than the case with no stenosis throughout the trachea. For 75% and 84%, there is a gradual increase in turbulence intensities up to approximately two diameters downstream of stenosis after which there is a smooth decay. However, it is interesting to note that the turbulence intensity for 84% constriction decays a little faster than that of 75%. For 91% constriction, we see a rapid growth in turbulence intensities for 15 l/min (Fig. 10.a) which peaks early at about half diameters downstream followed by a gradual fall. For 30 l/min, the kinetic energy continues to grow up to two diameters downstream before it starts decreasing which is due to the higher mass flow rate.

(13)

(a) (b)

Figure 11: Area-averaged pressure drop (p − pin) taken at every 5mm along the oral airway model. (a)

15 l/min (b) 30 l/min

(a)

Figure 12: Total pressure drop (pin−pout) for 15 l/min and 30 l/min

(14)

stenosis and explains why patients usually do not experience a major breathing impair-ment, or associated need for a stenting procedure (mechanical dilation of the airway), until the constriction is well above 50%.

(a) (b) (c)

(d) (e) (f)

Figure 13: Velocity profiles for a flow rate of 30 l/min with heliox as working fluid for 0% and 91% stenosis. (a,d) 2D streaklines. (b,e) Axial velocity contour (magnitudes in cm/s) and secondary velocity vector lines at A1-A2. (c,f) Axial velocity contour and secondary velocity vector lines at B1-B2

(15)

5 CONCLUSION

To the author’s knowledge, this is the first attempt to provide comprehensive character-ization of airflow in human oral airways with tracheal stenosis. For sedentary breathing (15 l/min), the flow remained laminar up to 75% stenosis. For normal breathing (30 l/min), even though the flow is turbulent with all stenotic constrictions, we noticed that the increase in pressure drops were only modest up to 75%, beyond which they rise sharply and are flow rate dependent. Compared to the idealized geometry, the evolution of la-ryngeal jet and the turbulence intensities were seen to be different for the case of normal breathing showing that the flow downstream of larynx is very sensitive to its shape. Since the airway geometry is completely case-dependent, it is important to consider realistic geometry in order to put CFD results to practical use. By using Heliox as an alternative breathing fluid, the effort of breathing is brought down by more than 50%.

6 ACKNOWLEDGEMENTS

Part of the research was funded by the VUB research council in the framework of a horizontal research activity and this funding is gratefully acknowledged. We thank Daniel Schuermans (Academic Hospital, Vrije Universiteit Brussel) for his assistance during the acquisition procedure of the CT scans.

REFERENCES

[1] J.R. Cebral and R.M. Summers. Tracheal and central bronchial aerodynamics us-ing virtual bronchoscopy and computational fluid dynamics. IEEE Transactions on Medical Imaging, 23, No.8, (August 2004).

[2] A.F. Heenan, E. Matida, A. Pollard and W.H. Finlay. Experimental measurements and computational modeling of the flow field in an idealized human oropharynx. Experiments in Fluids, 35, 70-84, (2003).

[3] G.M. Allen, B.P. Shortall, T. Gemci, T.E. Corcoran and N.A. Chigier. Computational simulations of airflow in an in vitro model of the pediatric upper airway. Journal of Biomechanical Engineering, 126, 604–613, (October 2004).

[4] G. Yu, Z. Zhang, R. Lessmann. Fluid flow and particle diffusion in the human upper respiratory system. Aerosol Science and Technology, 28, 146-158, (1998).

[5] T.B. Martonen, L. Quan, Z. Zhang and C.J. Musante. Flow simulation in the human upper respiratory track. Cell Biochemistry and Biophysics, 37, 27-36, (2002).

[6] C. Renotte, V. Bouffioux and F. Wilquem. Numerical 3D analysis of oscillatory flow in the time-varying laryngeal channel. Journal of Biomechanics, 33, 1637-1644, (2000). [7] I.M. Katz and T.B. Martonen. Flow patterns in the three-dimensional laryngeal

(16)

[8] M. Brouns, S. Verbanck and C. Lacor. Influence of glottic aperture on the tracheal flow. Journal of Biomechanics, (October 2005).

[9] T.E. Corcoran, N. Chigier. Characterization of the laryngeal jet using phase doppler interferometery. Journal of Aerosol Medicine, 13, 125-137, (2000).

[10] K. Wassermann, A. Koch, A. Warschokow and F. Mathen. Measuring in citu central airway resistance in patients with laryngotracheal stenosis. Laryngoscope, 109 1516-1520, (1999).

[11] V. Fasano, L. Raiteri, E. Bucchioni, S. Guerra, G. Cantarella, M.G. Massari, B.M. Cesana and L. Allegra. Increased frequency dependence of specific airway resistance in patients with laryngeal hemiplegia. European Respiratory Journal, 18 1003-1008, (2001).

[12] N. Hakimi. Preconditioning methods for time dependent Navier-Stokes equations. PhD Thesis, Dept. Fluid Mechanics, Vrije Universiteit Brussel, (1997).

[13] E. Turkel. Preconditioning methods for solving the incompressible and low-speed compressible equations. Journal of Computational Physics, 72, 277–298, (1987). [14] D. Choi and C.L. Merkle. Prediction of channel and boundary-layer flows with a

low-Reynolds number turbulence model. AIAA, 23, 1518–1524, (1985).

[15] A. Jameson, W. Schmidt and E. Turkel. Numerical simulation of Euler equations by finite volume methods using Runge-Kutta time stepping schemes. AIAA, 81-1259, (1981).

[16] C. Hirsch, C. Lacor, C. Rizzi, P. Eliasson, I. Lindblad and J. Hauser. A multi-block/multigrid code for the efficient solution of complex 3D Navier-Stokes flows. European Symposium on Aerodynamics for space vehicles, 415-420, (1991).

[17] Numeca’s Fine/Hexa user manual. version 2.1-a, (2005).

[18] C. Kleinstreuer and Z. Zhang. Laminar-to-turbulent fluid-particle flows in a human airway model. International Journal of Multiphase Flow, 29, 271-289, (2003). [19] Y.S. Cheng, Y. Zhou, B.T. Chen. Particle deposition in a cast of human oral airways.

Aerosol Science Technology, 31, 286-300, (1999).

[20] S.A. Ahmed and D.P. Giddens. Velocity measurements in steady flow through ax-isymmetric stenoses at moderate Reynolds number. Journal of Biomechanics, 16, 505-516, (1983).