## Abstract

Patients referred for treatment of tracheal stenosis typically are asymptomatic until critical narrowing of the airway occurs, which then requires immediate intervention. To understand how tracheal stenosis affects local pressure drops and explore how a dramatic increase in pressure drop could possibly be detected at an early stage, a computational fluid dynamics (CFD) study was undertaken. We assessed flow patterns and pressure drops over tracheal stenoses artificially inserted into a realistic three-dimensional upper airway model derived from multislice computed tomography images obtained in healthy men. Solving the Navier-Stokes equations (with a Yang-shih k-ε turbulence model) for different degrees of tracheal constriction located approximately one tracheal diameter below the glottis, the simulated pressure drop over the stenosis (ΔP) was seen to dramatically increase only when well over 70% of the tracheal lumen was obliterated. At 30 l/min, ΔP increased from 7 Pa for a 50% stenosis to, respectively, 46 and 235 Pa for 80% and 90% stenosis. The pressure-flow relationship in the entire upper airway model (between mouth and end of trachea) in the flow range 0–60 l/min showed a power law relationship with best-fit flow exponent of 1.77 in the absence of stenosis. The exponent became 1.92 and 2.00 in the case of 60% and 85% constriction, respectively. The present simulations confirm that the overall pressure drop at rest is only affected in case of severe constriction, and the simulated flow dependence of pressure drop suggests a means of detecting stenosis at a precritical stage.

- pressure drop

flow dynamics in the upper airways has gained a renewed interest in recent years, with the possibilities offered by computational fluid dynamics (CFD) to compute three-dimensional (3-D) particle deposition patterns for the benefit of drug therapy (21, 27). Nevertheless, the basic flow patterns of the carrier gas as it moves in from the mouth through a relatively complex bending oropharynx and into the trachea are still under study (2, 6, 10). It has been suggested that for flows as low as 3 l/min (21), turbulent flow patterns could develop because of local 3-D structural complexity. Different turbulence models have been applied to upper and tracheal airway structures, with varying degrees of geometrical simplification, to conclude that the flow regime is at least transitional between laminar and turbulent during resting breathing conditions, and turbulent during exercise (8, 16). For the present study, we use CFD simulations in a multislice computerized tomography (CT)-derived upper and tracheal airway structure to compute the effect of tracheal stenosis on the flow field. Our interest is less directed toward local flow patterns (which are crucial for local particle deposition) than toward overall pressure drops over the airway structure (which is a determinant of the work of breathing in stenosis patients).

Patients with tracheal airway stenosis often report a relatively sudden appearance of breathing impairment, which at the stage of admission to the clinic is observed when a loss of 75% or more of the airway lumen has occurred (20). Interventional bronchoscopic techniques such as laser or electrocautery resection with mechanical debulking of obstructive tissue, usually followed by airway stenting, then often have to be performed in a relatively urgent setting. There obviously is time for a progressive increase in tracheal stenosis before clinically significant breathing impairment is experienced by a patient, yet, when the constriction reaches a certain value, the pressure drop suddenly becomes critical and symptoms rapidly occur. With CFD simulations of upper airway flow dynamics in a realistic upper and tracheal airway structure, we aim to provide a better understanding of this clinical observation and to suggest a means of monitoring patients who are at risk for development of tracheal stenosis (e.g., patients with a history of longstanding or complicated intubation, posttracheostomy patients, lung cancer patients, or patients with a history of prior tracheal stenosis treated with laser or stenting).

A simple way of characterizing flow fields in a trachea or upper airway is to study its pressure-flow relationship. Wheatley at al. (25) measured in situ nasal and oral pressure drops in healthy men and related these to the corresponding flows to study the exponents of the power law that related both. When the power of flow had a value of 1.00, 1.75, and 2, these were assumed to represent laminar, turbulent, and orifice flow regimes, respectively. While the exact labeling of such categories may be debatable in strict fluid dynamic terms, the use of the power exponent to somehow quantify the flow regime at hand is interesting. Whether it can also provide a diagnostic tool that can be used in a clinical setting may be more complicated. In its most invasive experimental setting, this would require in situ pressure measurements beyond the stenosis, as was done by Wassermann et al. (24) in view of proposing a cut-off for stenotic resistance, beyond which surgery was indicated. A less invasive possibility is the use of frequency dependence of resistance, measured by body plethysmography, as was proposed by Fasano et al. (9) in patients with laryngeal hemiplegia.

The exploration of any of these possibilities also requires a detailed study of the flow fields in a representative airway structure with and without stenosis. The focus of the present study was on the so-called “weblike” stenosis, with a stenosis length of only 2 mm along the tracheal axis and typically occurring in the upper third of the trachea; this is a typical presentation of type I postintubation stenosis (19). For a subset of conditions, an elongated stenosis (length: 30 mm of tracheal axis) was also considered; this is a typical presentation of a “complex” or type II stenosis (19). With a particular attention to the relationship between pressure drop and flow over stenosis constrictions, we set out to quantify the exponents of the power laws at hand but also to determine the pertaining proportionality constants for easy estimation of actual pressure drops.

## MATERIALS AND METHODS

#### Model geometry.

Local ethics approval was obtained for a multislice CT imaging study in five otherwise healthy never-smoker male subjects who were referred for lung imaging due to possible low-concentration occupational asbestos exposure. After written informed consent was obtained, CT images were acquired using a Sensation 16 CT-scanner (Siemens, Erlangen, Germany) with the following settings: acquisition 16 × 0.75 mm, reconstruction slice thickness 1 mm, pitch 1.25, rotation time 0.5 s, care dose, reference mAS 100, 120 kV, medium filter B45f. Imaging was initiated above the nasal cavity and synchronized to start with the subject slowly inhaling to near total lung capacity where a breath hold was inserted to complete data acquisition just below the carina. In this way, the glottic opening was imaged with the subject being in a slow inhalation maneuver. From the scans obtained on the five subjects, one representative image of upper and tracheal airway anatomy was selected by the pulmonary radiologist and pneumologist.

The multislice CT images were imported, and a raw upper and tracheal airway 3-D geometry was constructed by triangulation (Amira, Mercury Computer Systems, Chelmsford, MA). From this raw geometry a smoothed geometry was derived, which included all features that were crucial to preserving the same flow field as in the raw geometry but that was better suited for artificially inserting stenoses of varying constriction and length. In a preliminary study where CFD simulations were used to guide particle imaging velocimetry experiments (3), flow fields in the raw and the smoothed realistic airway model had indeed shown excellent agreement in terms of maximum velocity, main recirculation zones, and location of the tracheal jet. The realistic (smoothed) model used here is shown in Fig. 1 with values for cross sections in Table 1, and its unstructured hexahedral mesh (Hexpress, Numeca, Brussels, Belgium) contains ∼650,000 cells with a local mesh refinement around the area of stenosis. A preliminary grid resolution study (using meshes with 433,000, 650,000, and 1,580,000 cells) had indicated that the velocity profile just downstream of the stenosis was not altered by a grid refinement above 650,000 cells and that the pressure drop between model inlet and outlet varied by <1% between simulations on the 650,000 and 1,580,000 meshes.

#### Computation methods.

The governing equations for the flow are the Navier-Stokes equations, and the available code (Fine/Hexa, NUMECA, Brussels, Belgium) is based on a compressible formulation and solves the Reynolds-averaged Navier-Stokes equations (RANS). The basic RANS equation integrated over a control volume V, with an elemental surface area d*S⃗*, is expressed as (1) where *U* is the conservative variable matrix containing density (ρ), the three components of mass flow (ρ*u*, ρ*v*, ρ*w*) and energy (ρ*E*); and *F⃗* and *Q⃗* are the inviscid and viscous fluxes, respectively. To close the system of RANS equations, a turbulence model is needed. In the present application a low Reynolds number Yang-shih k-ε turbulence model is used, the accuracy of which has been previously tested for the purpose of a similar application (2).

For detailed description of the computational methods, we refer to Brouns et al. (2). Briefly, because of the relatively low flow speeds in the present application, flow can be considered as incompressible, and solving such a problem with a compressible code leads to a stiff system with slow convergence. Hence, preconditioning is used. The preconditioning matrix employed here is a combination of those suggested by Choi and Merkle (5) and Turkel (22). A standard central scheme with Jameson type dissipation of second and fourth order is used for spatial discretization (15), and a fourth-order Runge-Kutta scheme for time integration. A full multigrid V-cycle strategy (11) with four grid levels was used for convergence acceleration. One typical simulation of a flow field in the realistic geometry of Fig. 1 took approximately 12–15 h to reach convergence (HP Integrity rx5670 4 × 1.3 GHz Itanium2, Hewlett-Packard, Palo Alto, CA).

## RESULTS

Figure 2 illustrates the flow fields for 30 l/min in the model geometries without stenosis (Fig. 2*A*) and with two degrees of a weblike stenosis (Fig. 2, *B* and *C*). The degree of constriction is quantified here as the obstructed cross-sectional area as a percentage of the nominal (unobstructed) airway cross-sectional area. The model *insets* in Fig. 2 show selected 3-D streamlines to illustrate the increasing complexity of flow field with increasing constriction. The main models in Fig. 2, *A–C*, show streaklines (i.e., the intersection of 3-D streamlines with a central plane), and the dark gray areas delimit areas where flows exceed 80% of peak velocity generated anywhere in the model. For instance, in the case of 0% and 50% stenosis, the highest velocities (∼5.5 m/s) are essentially in the vicinity of the jet generated by the glottic constriction, whereas in the case of 90% constriction, the peak velocities (∼22.5 m/s) are generated at the level of the stenosis and exceeding the velocities in the glottic area. As expected in the geometry without stenosis, the glottic area results in a so-called laryngeal jet, which impinges on the anterior side of the trachea with a huge flow separation zone forming at the posterior side of the trachea. The symmetry typical of a weblike stenosis basically reorients the jet along the tracheal axis, with flow separation occurring on either side of the jet in the case of a relatively mild stenosis (e.g., 50% constriction; Fig. 2*B*). In the case of a more severe constriction (e.g., 90% constriction; Fig. 2*C*), the laryngeal jet is dwarfed by the jet developed at the site of stenosis.

In Fig. 3, pressure averaged over the cross sections along the upper airway axis is plotted as a function of distance from the model inlet for a flow of 30 l/min. The effect of the glottic narrowing (located at 152 mm from model inlet) is seen to be of the same order as that of the 50% constriction, 34 mm further downstream. By contrast, a stenosis narrowing to 75% constriction (i.e., leaving 25% of the tracheal lumen unobstructed) more or less doubles the pressure drop across the stenosis. Extra 10% and 15% constrictions produce dramatic increases in pressure drop. For the greatest stenosis constriction (90%), the pressure obtained in longer stenoses with a central channel of 10, 20, and 30 mm length is also represented. These long stenoses show a modest difference in pressure profile with a slightly smaller magnitude of total pressure drop than the weblike stenosis of comparable constriction (90%). The difference in pressure drop between the 10- and 30-mm stenosis (26 Pa) is approximately twice the pressure drop between the 10- and 20-mm stenosis (12 Pa).

When expressing the simulated pressure drops over the stenosis as a function of constriction (Fig. 4*A*), a dramatic ΔP increase with constriction is observed well beyond 70% constriction, in a characteristic pattern that is similar for 15 and 30 l/min inlet flows. When considering heliox instead of air at a flow rate of 30 l/min (Fig. 4*A*), this resulted in a reduction of pressure drop to approximately one-third of its value with air at 30 l/min for all constrictions in Fig. 4*A*, consistent with the density of heliox being one-third that of air.

From these CFD-simulated curves, we attempted to derive an empirical formula that could serve as a rule of thumb in a clinical situation where inlet flow and stenosis cross section can be reliably estimated. For this purpose, we first consider the generalized Bernoulli equation (excluding the effect of gravity), which takes into account cross-sectional changes (as if they were occurring in a straight tube) and a sum of so-called losses ΔP _{i}^{loss} that account for frictional losses (e.g., induced by the overall geometry) and local losses owing to particular geometrical features (e.g., an orifice at the level of the stenosis) (2) where P is pressure, *indexes 1, 2,* and *3* refer to locations at, respectively, model entry, stenosis, and model exit; ρ is gas density (1.2 kg/m^{3} for air; 0.4 kg/m^{3} for heliox), *V* is local velocity, and *S* is local cross section. If we consider that in our particular case, the term ΔP_{i}^{loss} refers to a simple sum of the pressure loss generated by the overall upper airway geometry (ΔP_{[no stenosis]}) and the pressure loss generated by the stenosis (ΔP_{[stenosis]}), the latter term can be isolated to estimate the pressure drop generated by the stenosis. In engineering applications, ΔP^{loss} (and applied here to ΔP_{[stenosis]}) is usually expressed as (3) where Q is bulk flow and *K* a constant that is empirically determined for each particular geometrical feature (in this case, location 2 at the stenosis); *K* = 0 corresponds to zero loss, and a typical *K* value for airflow in a 90° miter bend is 1.1, but *K* can increase to more than 100 for a partially opened butterfly valve (26).

By fitting the simple *Eq. 3* to the CFD-simulated ΔP over the stenosis represented in Fig. 4*A*, *K* was found to vary between 0.4 and 1.5 for 15 and 30 l/min flows and for constrictions between 50 and 90% (Fig. 4*B*). For the sake of simplicity, with the aim to turn *Eq. 3* into one rule of thumb that is valid for all conditions within the above range, we chose one single value for *K* that provided a good fit for the 15 and 30 l/min curves in Fig. 4*A*. The solid line (•) and dotted line (○) in Fig. 4*A* are pressure drops for, respectively, 15 and 30 l/min, obtained from *Eq. 3* with *K* = 1.2, plotted alongside the 15 and 30 l/min CFD-simulated data points for air and heliox.

The general form of *Eq. 3*, typically applied in engineering problems, considers a dependence on bulk flow with a power law and exponent 2, which is a good approximation for a rule of thumb designed to roughly estimate a local pressure drop generated by a specific feature (stenosis) as a function of bulk flow. In the present context, it was useful to also perform a more detailed inspection of pressure drop dependence on flow, considering the entire upper and tracheal airway model of Fig. 1. Figure 5 shows the pressure drops obtained between model inlet and outlet for different flows up to 60 l/min, in the case of no stenosis and of 60% and 85% constriction. On the represented data points we fitted a power law (ΔP = *aV̇*^{b}, where *V̇* is flow rate) to obtain an exponent *b* to characterize pressure dependence on flow. In the presence of 60% and 85% constriction, the best power fit exponent *b* was 1.92 and 2.00 (*R*^{2} = 1.00 for both). In the absence of stenosis, the best fit *b* was 1.77 (*R*^{2} = 0.999). When using only data points below 30 l/min, *b* was similar (1.72; *R*^{2} = 0.999), but a linear regression forced through zero in this flow range up to 30 l/min led to a worse fit with *R*^{2} = 0.89.

## DISCUSSION

In the present study, CFD-simulated pressure drops over tracheal weblike stenoses of different constriction in a realistic upper and tracheal airway geometry roughly follow a dependence on the squared ratio of bulk flow to stenosis lumen cross section. This relationship predicts a dramatic and relatively sudden increase in breathing impairment experienced by patients once a relatively severe degree of stenosis is reached, in agreement with what is observed in clinical practice. This characteristic dependence of pressure drop on increasing degrees of weblike stenosis allowed us to derive a simple equation to actually quantify pressure drops generated by different degrees of stenosis (*Eq. 3*), provided that breathing flow and stenosis cross section can be reliably estimated. The coefficient *K* in *Eq. 3* is expected to partly depend on the development of the velocity profile upstream from the stenosis and on the exact geometry of the constriction, explaining slightly varying *K* values between 15 and 30 l/min, and for different degrees of constriction (Fig. 4*B*). However, the agreement between the CFD simulations and *Eq. 3* by use of a single *K* value (*K* = 1.2; Fig. 4*A*) indicates that this rule of thumb can be adequately used for most clinically relevant situations (slow to quiet breathing).

For the purpose of decision making on weblike stenosis treatment, pressure drops owing to stenosis can be roughly estimated by means of *Eq. 3* with *K* = 1.2, given that a noninvasive measure of the stenosis cross section is available, for instance by an acoustic reflection technique (12). The resulting resistance at a given flow can then be computed and related in terms of abnormality with respect to a reference value of upper airway resistance. In Wassermann et al. (24), an inspiratory resistance sevenfold greater than that measured in normal control subjects was used as a cut-off to discriminate between stenosis patients who needed surgical treatment or not; for this purpose Wasserman et al. (24) measured tracheal pressure drops in situ. In fact, our 15 and 30 l/min simulations show that a sevenfold increase in pressure drop would already correspond to 85–90% constriction (Fig. 4*A*). Irrespective of the preferred choice for a cutoff, the present study suggests that a similar approach could be used but by estimating pressure drops from *Eq. 3*, where only stenosis cross section needs to be measured, which may be of considerable practical advantage with respect to actually having to measure the pressure drops in situ (24).

A more detailed inspection of the flow dependence of pressure drop over the entire model provides another perspective for clinical decisionmaking with respect to tracheal stenoses. Indeed, the dependence of pressure drop on flow, with its specific power law exponent depending on whether a stenosis is present or not (Fig. 5), suggests an attractive noninvasive way to probe the trachea for presence of considerable stenosis. Wheatley et al. (25) have determined the power law exponent to quantify pressure-flow relationships they were measuring in human nasal and oral passages, based on the rationale that power exponents 1.00, 1.75, and 2 would, respectively, correspond to laminar, turbulent, and orifice flow regimes. Recording oral and nasal pressure-flow power relationships, they obtained best-fit exponent values around 1.75, with a consistently greater exponent during inspiration (1.83 ± 0.04) vs. expiration (1.66 ± 0.07), which was thought to signify greater level of turbulence during inspiration. Surely, the high quality of fit in Wheatley et al. (25) was partly due to the fact that pressure and flow were measured in situ, over very well-defined oral or nasal pathway structures. Yet, even if it is not possible to directly access pressure in situ in the case of tracheal stenosis, one could imagine that tracheal contribution to a noninvasive measure of upper resistance would suffice to identify a characteristic flow dependence of resistance. An increase of the best-fit power law exponent between 1.75 and 2 during a given patient's followup could then signal a progression toward stenosis that impairs breathing at rest, especially because the exponent is seen to increase to 1.92 even in the case of 60% constriction, which in Fig. 4*A* corresponds only to a preliminary stage of ΔP increase. In any case, Fig. 5 clearly illustrates how resistance measurements at higher flows could be more useful as an indicator of flow limitation due to stenosis (see, for instance, a 50% increase in resistance between no stenosis and 60% stenosis when measured at 1 l/s).

Little is known about pulmonary function measurements in stenosis patients that could actually lead to a better management of stenosis treatment. Because of the rapid onset and progression of symptoms in most patients with tracheal stenosis, pulmonary function testing including spirometry and flow-volume loop analysis can rarely be obtained in patients referred in extremis. Diagnosis severity assessment and therapeutic decision in tracheal stenosis are therefore based on history, physical examination, imaging studies, and flexible bronchoscopy. The role of pretreatment pulmonary function testing outside of research protocols is unclear (7). Significant improvement in forced vital capacity, forced expiratory flow in 1 s, and in peak expiratory flows, forced inspiratory volumes, and airway resistance has nonetheless been demonstrated after airway stenting (23). There is no standardized approach in the followup of patients treated for, or at risk for, tracheal stenosis. Timing of followup visits, including history and clinical examination, imaging, pulmonary function testing, and/or flexible bronchoscopic control, should be individualized; the role of surveillance bronchoscopy has not been established (17). The present study suggests that flow dependence of upper airway resistance measurement (Fig. 5) could provide a valuable diagnostic tool. In addition, the heliox simulations (Fig. 4*A*) also lend support to the experimental observation of a temporary relief of breathing impairment during heliox administration as has been suggested for use with postintubation patients (13).

The glottic constriction of the tracheal model used here was 125 mm^{2}, which corresponds to the average 126-mm^{2} peak value of the glottic area obtained by Brancatisano et al. (1) at midinspiration. The difference in pressure drop obtained here for 30 l/min (15 Pa) compared with those obtained on simplified models (2, 10) (∼30 Pa across both studies) can be almost entirely accounted for by the difference in glottic cross section. Given that in these previous models, cross sections ranged from 0.9 to 1 cm^{2}, and considering that simulated pressure drop was seen to decrease approximately as a square of increasing cross section (2), the squared ratio (0.95/1.25) = 0.58 accounts for most of the difference in pressure drop between the present model without stenosis and previous unobstructed models. The remaining difference can be probably accounted for by the abruptness of cross-sectional changes in the oropharyngeal structure proximal to the glottis. In any case, the present study clearly shows that the pressure drop over the glottis that actually corresponds to a 40% constriction (i.e., 125 mm^{2}/310 mm^{2}; see Table 1), is negligible with respect to those induced by the constrictions greater than 70%, which impair breathing. In fact, Jaeger and Matthys (14) used an equation of the form of *Eq. 3* to study the effect of laryngeal constriction on the relationship between pressure drop and flow. Using the plot of the experimentally obtained discharge coefficients (corresponding to 1/*K* in *Eq. 3*) vs. Reynolds number [as illustrated in Fig. 11 of Chang (4)], and assuming that the same relationship would hold for stenosis constrictions, we computed the estimated ΔP for the cases considered in Fig. 4*A*. For 50% constriction, ΔP values were similar to those obtained with CFD in our model, whereas for 80% and 90% constrictions, the CFD-computed ΔP was underestimated by, respectively, 30% and 50% when using the discharge coefficient proposed for the larynx.

While the focus of this work was on weblike stenosis, we did investigate the effect that an elongated stenosis could have on pressure drop in the case of the most severe constriction. In fact, the pressure drop with the three long stenoses is slightly smaller compared with that obtained for a weblike stenosis of similar constriction (90%). This is due to the fact that in the weblike stenosis, the constriction is immediately followed by the expansion, while in the long stenoses the flow is able to develop in the small channel between constriction and expansion. The dependence of pressure drop on length of the long stenosis can be explained on the basis of the frictional losses that come into play along the channel lengths. In turbulent regime, the pressure drop due to friction of a fluid in a horizontal pipe of length *L* and diameter *d* can be approximated by 0.158*L*ρ^{0.75}μ^{0.25}*d*^{−1.25}*V*^{1.75} (18, 26), where ρ is density of the fluid, μ is dynamic viscosity, and *V* is velocity in the tube; the coefficient 0.158 was taken from White et al. (26), which was considered representative of a rough-walled tube. In the case of 90% stenosis, flow rate of 30 l/min and a stenosis length of 10, 20, and 30 mm, respective frictional losses are then estimated at 8.6, 17.1, and 25.7 Pa (i.e., corresponding to the multiplicative factor on length). Considering the difference in inlet conditions to the stenosis of our model vs. that of a straight tube of diameter *d* and length *L* to which this equation applies, the CFD-computed pressure drops between 10- and 30-mm stenosis (26 Pa) are more or less double that between 10 mm and 20 mm (12 Pa), which is consistent with the above-predicted frictional losses.

In conclusion, by applying CFD in a realistic model of tracheal stenosis, we simulated pressure drops during normal breathing that are consistent with the clinically observed sudden appearance of breathing impairment at a relatively advanced stage of tracheal obliteration. A rule of thumb is derived from which pressure drops over the stenosis can be estimated on the basis of breathing flow and stenosis cross section. In addition, the CFD simulations enabled us to suggest that the best-fit exponent in the power law that relates upper airway resistance to flow could indeed be used as a diagnostic tool in the noninvasive monitoring of stenosis patients.

## GRANTS

This study was funded by the Fund for Scientific Research-Flanders. A part of this research was funded by the Vrije Universiteit Brussel Research Council in the framework of a horizontal research activity and this funding is gratefully acknowledged.

## Acknowledgments

We thank Daniel Schuermans (Academic Hospital, Vrije Universiteit Brussel) for assistance during the acquisition procedure of CT scans.

## Footnotes

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

- Copyright © 2007 the American Physiological Society