The effects of curvature and constriction on airflow and energy loss in pathological tracheas

This paper considers factors that play a significant role in determining inspiratory pressure and energy losses in the human trachea. Previous characterisations of pathological geometry changes have focussed on relating airway constriction and subsequent pressure loss, however many pathologies that affect the trachea cause deviation, increased curvature, constriction or a combination of these. This study investigates the effects of these measures on tracheal flow mechanics, using the compressive goitre (a thyroid gland enlargement) as an example. Computational fluid dynamics simulations were performed in airways affected by goitres (with differing geometric consequences) and a normal geometry for comparison. Realistic airways, derived from medical images, were used because idealised geometries often oversimplify the complex anatomy of the larynx and its effects on the flow. Two mechanisms, distinct from stenosis, were found to strongly affect airflow energy dissipation in the pathological tracheas. The jet emanating from the glottis displayed different impingement and breakdown patterns in pathological geometries and increased loss was associated with curvature.


Introduction
Many pathologies affect the flow of both air and blood in the human body by changing the shape of the relevant conduit. This change of shape has often been characterised by a constriction to the conduit (e.g. degree of stenosis) and related to the pressure drop. Several studies have hypothesised that the degree of constriction can be used as a predictor for the physiological impact of the pathology. For example, in the airways, Brouns et al. (2007) considered the degree of tracheal stenosis in a straight tracheal geometry, whilst in haemodynamics Schrauwen et al. (2015) investigated pressure losses in straightened atherosclerotic coronary arteries. However, the constriction of the conduit may be only one of several geometric changes caused by the pathology. Typically, conduits in the human body exhibit some degree of curvature and certain pathologies increase curvature. These additional geometric factors have previously been overlooked in the airway literature. Brouns et al. (2007) presented a means to predict increased tracheal pressure loss from the degree of constriction. Their study introduced an index to calculate pressure loss from stenosis degree, however the analysis did not include the influence of curvature (or other geometric changes) thus the index falls short of a complete geometric characterisation. The importance of different geometrical effects was highlighted by Mylavarapu et al. (2013) who investigated the resistance in a trachea with a subglottal stenosis. They focused on changes to this resistance in four possible surgical procedures. All four virtual surgery cases increased the cross sectional area of the airway, yet two cases demonstrated a rise rather than the expected fall in airway resistance, demonstrating that constriction is not the only factor in determining tracheal resistance. This study highlights the need for a more in depth understanding of the relationship between airway anatomy and the impact this has for the subject in terms of resistance to airflow and breathing effort, particularly in terms of how surgical procedures may modify this anatomy.
Analysis of tracheal flow mechanics using CFD has seen significant attention in healthy geometries (Lin et al., 2007;Choi et al., 2009;Ma and Lutchen, 2006;Comerford et al., 2013). These studies have observed a number of fluid mechanics phenomena as the flow enters the trachea through the constriction of the glottis, these include: breakdown and impingement of the jet; turbulence; and secondary flow vortices due to the curvature of the vessel. However the trachea can be affected by a number of pathologies that further influence the geometry causing constriction, deviation or a combination of both, thus altering the flow mechanics and breathing effort.
A number of studies have investigated the influence of different pathologies on tracheal flow mechanics (Brouns et al., 2007;Mihȃescu et al., 2008;Malvè et al., 2011;Mimouni-Benabu et al., 2012;Hamilton et al., 2015). As already mentioned, Brouns et al. (2007) simulated the influence of tracheal stenosis in an idealised tracheal geometry. The stenosis they considered was axissymmetric and positioned in the sub-glottal region. Their results indicated that for an idealised stenosis, the relationship between pressure and flow (a power law) increases from normal only when severely constricted (60% and 85% luminal area reduction).
Further studies on flow in stenosed airways have been performed by Mihȃescu et al. (2009), Mimouni-Benabu et al. (2012 and Vial et al. (2005). These studies further highlight the need for increased understanding of how airway anatomy affects the work of breathing and how this can affect clinical and surgical planning.
The degree of involvement of the upper airways in defining flow needs to be considered. Lin et al. (2007), Xi et al. (2008), Choi et al. (2009) andPollard et al. (2012) have demonstrated that including some of the upper airway geometry is necessary for an accurate description of the fluid mechanics in the trachea. This is primarily because the glottal jet -which forms due to the glottal constriction -influences turbulent flow structures in the trachea. As to how much of the upper airways must be considered, Saksono et al. (2011) reported that the whole nasal cavity must be included, however Choi et al. (2009) reported that simulations truncated at the mid-laryngopharynx level provide a good description of tracheal flow mechanics (when compared to simulations of the whole upper airways), with similar results found with truncation immediately above the glottis. Therefore, inclusion of the most constricted section upstream of the trachea is deemed sufficient -in this study, this is assumed to be the glottis.
Previously, LES approaches have been presented for modelling airflow in healthy subjects (Lin et al., 2007;Choi et al., 2009;Comerford et al., 2013). Mihȃescu et al. (2009) compared flow calculated using steady Reynolds-averaged Navier-Stokes (RANS) calculations and unsteady LES in a paediatric trachea with subglottal stenosis, concluding that LES is the preferred tool to capture the flow features associated with the stenosis such as large radial velocity gradients. In the study reported here, the accuracy of LES is compared with that of very high resolution simulations devoid of turbulence modelling assumptions. CFD simulations have also been found to closely match experimental results in an idealised healthy extra-thoracic airway (Heenan et al., 2003;Ball et al., 2008). Some differences were observed, but can potentially be attributed to limits on mesh resolution. Whilst further previous studies have investigated air flow in stenotic tracheas, those not including the glottis are not considered here. Calmet et al. (2016) and Bates et al. (2015a) both analysed the flow regime in highly resolved simulations, concluding that flow in the complex geometries of the airways at these flow rates is likely to be non-laminar with mild turbulence observed.
Highlighting the increase in respiratory effort through properties other than constriction has clinical significance. The American Thyroid Academy (Stang et al., 2012) propose basing the surgical decision making on the ratio of 1D measurements (the airway diameter in a CT slice) of the tracheal anatomy, a measure which again does not consider the influence of other geometric factors. In this study we investigate the relative importance of these previously overlooked geometric measures through one specific pathology, the compressive goitre.
Combining geometric measures with CFD as a clinical management tool for pathological changes has been demonstrated by Sagittal and coronal views of tracheal geometries. A was described as normal, B and C show both curvature and constriction, whilst D and E show increased curvature from the normal case, but without constriction. Each geometry extends above the glottal region shown here to ensure that natural inflow conditions have developed across the area of interest. The dashed white line represents the location of the first tracheal ring. Hamilton et al. (2015), who used area and pressure loss as measures for functional analysis of the development of a transplanted trachea.
A goitre is an enlargement of the thyroid gland, which in certain cases can extrinsically compress the trachea, causing it to displace or deform. This can lead to symptoms such as difficulty in swallowing, coughing and shortness of breath or an inability to breathe at higher flow rates. The decision to operate to remove a goitre is based on a number of patient parameters, including analysis of CT images. Additional diagnostic complications arise when patients suffer also from lung pathology, masking the influence of the compressive goitre.
The primary purpose of the present paper is to demonstrate how CFD simulations of tracheal flow can enhance understanding of the relationship between pathological tracheal anatomy and losses of pressure and energy. This relationship can provide useful insights for surgical planning. Given the limited number of cases considered, necessitated by the cost of detailed LES simulations, the focus is not to examine the broad range of possible geometries, but to facilitate a better understanding of how geometry affects the airflow and the consequences in terms of factors that must be considered in any clinically relevant metric.

Selection of geometry
To provide an in-depth analysis of geometric factors that are important for breathing mechanics five tracheal geometries were studied: one normal and four pathological cases (affected by compressive goitre). The normal case provides a reference, similar to the studies of a number of healthy tracheal geometries in the literature (Lin et al., 2007;Choi et al., 2009;Ma and Lutchen, 2006;Longest and Vinchurkar, 2007) to compare with the pathological cases.
The five geometries which were studied are shown in Fig. 1. All geometries are imaged at high lung volume, nominally total lung capacity, sustained as long as the subject was able. Images were then chosen with abducted vocal cords.
Geometry A was deemed normal by a consultant radiologist and is used for comparison both with the pathological cases and with other normal cases in the literature. (Subsequent analysis of the flow within this case exhibits similar values as other normal tracheas described in the literature across a range of key flow parameters, as described in Section 3.1, thus validating its selection as a reference case.) Geometries B to E are all deemed to display a compressive goitre, however they were chosen as they fall into two distinct classes. Geometries B and C exhibit both constriction and curvature above the normal range. Geometries D and E are not constricted, as defined later in Section 2.3. However, both exhibit elevated curvature.
These geometries were chosen as they allow comparison between the effects of curvature, which has not previously been assessed in tracheal anatomies both with and without the effect of constriction, a factor which is known (Brouns et al., 2007) to increase losses at high degrees of stenosis.
Realistic geometries, derived from medical CT images, were chosen rather than idealised geometries due to the complex nature of the anatomy of the larynx and the many effects it has on the airflow within. Several previous studies (Lin et al., 2007;Choi et al., 2009;Ma and Lutchen, 2006;Comerford et al., 2013) have studied the flow through realistic healthy laryngeal geometries and have identified the dominant flow features; the angle and point of separation of the laryngeal jet, and the location and angle of its impingement on the tracheal wall. The effect of each of these significantly changes the flow patterns within the airway (Bates et al., 2015b). Idealisation of the airway anatomy may easily neglect the complex relationship between these flow features and the anatomy (Banko et al., 2015), hence realistic geometries are deemed more suitable.

Data acquisition
Anonymised CT scans of the neck and chest were obtained retrospectively from an image database. Ethical approval to use the data from the local joint research compliance office was waived and research and development approval was granted. Only scans containing the necessary anatomy and reported by a consultant radiologist as displaying a compressive goitre distorting tracheal anatomy from the first tracheal ring to the carina were selected. Clinical opinion was that the glottis was not affected by the goitre in any of the cases. Furthermore, only scans with a slice thickness less than 2.5 mm were chosen, to aid accurate surface extraction. All scans were segmented by a trainee ear, nose and throat (ENT) surgeon using Mimics 17 (Materialise NV, Leuven, Belgium). The segmentation process involves marking the voxels of the CT image set which represent the airway, a threshold value of −500 Hounsfield units was used to delineate the air-tissue boundary, with manual correction where necessary to generate a realistic physiological airway surface. This procedure produced 3D surfaces of the geometry which were then smoothed using Taubin's (Taubin, 1995) algorithm in VMTK 1.2 (Vascular Modelling Toolkit, http:// www.vmtk.org/) to eliminate the stepped nature of the voxel surface.

Geometric analysis
Centrelines were produced from the surface mesh of each geometry using the software VMTK, following the methods described by Piccinelli et al. (2009). The centrelines consisted of a series of equispaced points along the trachea. Using an in-house MATLAB 2014b (Mathworks, Natick MA, USA) code, cross-sectional planes were generated normal to the centreline and bounded by the geometries' surfaces. The areas and perimeters of these planes were then measured.
To understand geometric features in detail we utilise the Frenet reference frame, which carries an intrinsic description of the local geometry (Piccinelli et al., 2009). In this description the centreline is parameterised by a variable s that describes the arc length. At discrete positions along this centreline a position vector, c(s), is defined relative to the origin. From this position vector the following unit tangent (T(s), points along the centreline) and normal (N(s), points towards the centre of curvature) vectors can be defined: where Ä = c (s) is the curvature. The cross-sectional areas of the planes normal to the centreline provide the following metrics: A MEAN , the mean cross sectional area of the trachea from the first tracheal ring to one diameter above the carina, A GLOTTIS , the glottal area, the ratio between this value and the mean area of the trachea serves as an indication of the strength of the glottal jet and AR, the area reduction, representing the severity of the constriction and is defined by where A MIN is the minimum cross-sectional area, A 1ST RING is the cross-sectional area of the first tracheal ring and L TRACHEA is the length of the trachea centreline, measured from the 1st tracheal ring to the carina. Table 1 gives these geometric quantities for each case. There is naturally a degree of variation in the geometry of normal tracheas. Those with goitres exhibit further variation due to a pathological extrinsic compression and displacement resulting in abnormal curvature of the trachea.
In order to assess the extent of compression, current guidelines for clinical assessment make use of tracheal constriction deduced via measurement of airway diameter on an image plane. Specifically, the American Thyroid Association guidelines (Stang et al., 2012) suggest a 35% or larger reduction in diameter from the first tracheal ring to the narrowest point to indicate a pathological anatomy.
This study uses area reduction (AR) rather than diameter reduction to determine the degree to which each case is constricted, as cross-sectional area is a more thorough measure of the space available to the airflow than a 1D diameter. Constrictions due to goitres are not uniform in all directions, therefore we have used the value of 35% proposed by Stang et al. (2012) to divide our pathological geometries into two classes: those tracheas exhibiting variation in cross-sectional area within the normal range and those which are constricted by more than 35% of the first tracheal ring. Whilst this measure is a ratio of diameters, the non-axisymmetric nature of the constriction in goitre subjects means that the area reduction is likely to scale linearly with it. The classification of constriction is not particularly sensitive to this value of stenosis percentage, any value between 30% and 50% would have divided the cases into the same groups, as shown in Table 1. The table also shows that all four pathological cases have a mean curvature of at least double the value of the normal case.

Mesh generation and simulation
Volume meshes were generated using Star-CCM+ 8.04.007-r8 (CD-adapco, Melville, NY, USA) and the unsteady, incompressible Navier-Stokes equations were solved using the same software. Polyhedral cells were used throughout the bulk of the flow, as they can capture flow features at lower computational expense than an equivalent tetrahedral mesh (Peric, 2004). Prism layers were used to line the walls, depending on the mesh resolution. Full details of the solution strategy and validation can be found in Bates et al. (2016). All geometries were extruded at both the inlet (above the glottis) and at the outlets (below the carina) to ensure well defined boundary conditions. Flow measurements were only taken below the constriction at the glottis, to ensure that natural inflow characteristics had developed (Choi et al., 2009).
Quasi-steady simulations, that is unsteady simulations with constant boundary conditions, were performed as the best compromise between computationally expensive unsteady simulations and faster, steady-state simulations, which neglect much of the unsteady dynamics observed within tracheal airflow, such as fluctuations of the glottal jet. These quasi-steady simulations were deemed sufficient to determine the effect of geometric changes on inhalation as Bates et al. (2015a) showed that for 90% of the duration of the inhalation they considered, flow patterns remained constant, with velocity magnitudes scaling with the driving flow rate. Therefore, the majority of the duration and volume of an inhalation will be contained within this phase. Inhalation was chosen as the CT images were captured during inhalation, considered to be the limiting breathing phase for these subjects. The geometry of the large airway differs significantly between inhalation and exhalation (Schwab et al., 1993), hence exhalation should not be analysed for this dataset.
Using quasi-steady rather than steady simulations allows realistic temporal dynamics and fluctuations to be considered. Initially, the air is at rest throughout and flow is allowed to develop until start-up transients have attenuated, as described in Bates et al. (2016). The simulations are then continued until the temporal mean values have converged.
The boundary condition at the inlet was a uniformly distributed velocity profile. Whilst non-physiological, this boundary is applied far enough upstream of the region of interest to allow realistic flow to have developed through the glottis and the region downstream where flow is analysed (Choi et al., 2009). The flow at the distal end of the trachea, the outlet, is held at constant pressure, but the flow profile here is allowed to develop naturally. Flow rates of 15, 30 or 60 l min −1 were studied as they represent the peak instantaneous flow rates which may be achieved during restful breathing, manoeuvres such as sniffing, through to moderate exercise (Wasserman et al., 1973;Comer et al., 2001;Rennie et al., 2011).
The downstream boundary was simply held at a constant reference pressure. This study only considers pressure losses between sections of interest, so the absolute pressure value has no influence on the results. Again, whilst this boundary condition is non-physiological, it has no influence on flow in the region of interest (terminating one tracheal diameter above the bifurcation at the carina), the velocity and pressure fields are unaffected by the peripheral flow distribution (Tadjfar and Smith, 2004;Smith et al., 2003;Shah et al., 2012).
Using pressure driven boundary conditions rather than an applied flow rate is not possible as the pressure differential between the ends of the trachea is not easily measured in vivo, whereas the instantaneous breathing flow rate may be recorded, as it is in spirometry. Bates et al. (2015a) have repeated simulations with equivalent pressure and flow driven simulations and found no difference in the flow fields or quantities of interest.
The no slip condition was applied at all walls. These simulations took approximately 3500 CPU hours (per simulation) to solve, or three days on a 48 core cluster, whilst this may be too high for clinical use, these simulations were performed at relatively high spatiotemporal resolution in order to capture the majority of the fluid dynamics of interest.

Calculation of metrics
To determine the pressure drop along the trachea the mean total pressure was calculated on two cross sectional planes. The first plane was located at the level of the first tracheal ring, whilst the second was one diameter upstream of the apex of the bifurcation at the carina. The choice of suitable end point for measurement of tracheal loss is motivated by the theoretical and computational studies of Tadjfar and Smith (2004) and Smith et al. (2003). Whilst they found that pressure varied strongly in flow approaching a bifurcation or a side branch, the region of upstream influence is effectively limited to one diameter.
The tracheal resistance is then found from dividing the tracheal pressure loss by the flow rate.
Analysis of the influence of geometry on complex three dimensional flows within curved, non-symmetric ducts can be aided by reducing the data to a one dimensional (1D) measure along the length of the trachea. The distribution of flow in both stenotic (Varghese et al., 2007) and curved (Berger et al., 1983) geometries has been shown to be non-uniform. Furthermore, flow in stenosed tracheas is dominated by jets (Choi et al., 2009) and regions of recirculation (Mihȃescu et al., 2009). Each of these features reduces the proportion of the lumen's cross sectional area which carries the main through-flow of air. In order to quantify the cross-sectional area effectively used by the bulk of the flow, two measures of utilised area (UA) are proposed; UA MAX and UA MEAN . The area UA MAX is the section of the cross-sectional plane where the local plane-normal velocity, U ·n is above half the maximum U ·n on that cross-section (n being the normal unit vector to the plane). The parameter UA MEAN is the area where U ·n is greater than the mean value for that plane. Both UA MAX and UA MEAN are then expressed as ratios to the total available area. In Fig. 2 a sample airway lumen is shown with the region UA, corresponding to the part favoured by the flow marked.
The metric UA MAX quantifies the degree of flow non-uniformity, which is an inevitable consequence of curvature and constriction. As flow through an asymmetric constriction is associated with jet impingement and subsequent breakdown, these two effects redistribute flow in the azimuthal direction: UA MAX will be reduced in these locations. Distal to the breakdown region, the redistribution of kinetic energy leads to a flattening of the velocity profiles. UA MAX increases as the high velocity redistributes.
The energy flux, E through each cross-sectional plane was also calculated, where U is the velocity vector,n is the normal unit vector of the plane, which has the surface S. P s is the static pressure on the plane and is the density of air. The change in energy flux along the length of the trachea equals the power required to drive the flow through the geometry.

Limitations
This study considers quasi-steady simulations of inhalation. The quasi-steady assumption means the formation of flow patterns and their decay are not considered (Bates et al., 2015a), exhalation is also not considered due to incomplete information of the inflow condition from the lower airways. Also, the geometries were derived from medical images taken during inhalation and the airway geometry differs in each phase of breathing (Schwab et al., 1993). The walls are considered to be rigid in this study due to the lack of dynamic imaging of goitres. All scans were performed with the subject in the supine position during inhalation (whilst the subjects were awake), variations to the anatomy in other positions are not considered.

Overall assessment
A basic assessment criterion for the effort of breathing is the pressure loss; that is the pressure difference between the ends of the trachea that must be generated to drive a given flow rate. The current assumption, in wide use clinically, is that constriction is the dominant cause of elevated pressure loss in pathological cases. Higher pressure losses and resistances imply a higher work of breathing.
Comparing the work required to drive flow through the trachea to the total work of breathing, may indicate at what severity tracheal pathologies start to cause symptoms such as breathlessness.
The pressure loss in the five geometries above (which exhibit a range of constrictions) was calculated at a range of flow rates and analysed to determine the underlying functional relationship.
The pressure-flow relationship for each tracheal geometry is plotted in Fig. 3. The total pressure loss is the difference between the surface average of the temporal mean of total pressure on planes across the glottis and above the carina. The use of a logarithmic plot shows a simple power law relationship provides an excellent fit to the data within the flow range of interest and allows the corresponding exponent to be readily deduced. Table 2 gives the exponent for each geometry. Brouns et al. (2007) found that an idealised oral geometry exhibited a pressureflow exponent of 1.77, the same value as the healthy geometry considered here (A). However, Brouns et al. (2007) reported that the exponent increased with degree of constriction, whereas in this case, geometries B, C and E all have similar exponents, despite exhibiting different abnormalities: 55%, 60% and 17% stenoses respectively. Whilst all pathological cases have a higher exponent than that found in the normal case, the greatest value is found in the unconstricted (stenosis <35%), but highly curved geometry D. Both curved, but unconstricted geometries (D and E) have high exponents with a low resistance.
The constricted cases (B and C) have relatively high losses at low flow rates, but resistance increases more gradually as the flow rate rises. Scaling the pressure values by either the glottal area or the square of the glottal area does not influence pressure-flow characteristics, indicating that this behaviour is not merely a function of the constriction at the glottis. Table 2 also shows the resistance values from the glottis to the carina for each geometry. These values are calculated by dividing the pressure loss by the flow rates shown in Fig. 3. Both McRae et al. (1995) and Smitheran and Hixon (1981) report resting breathing rate resistances in the region of 0.006 Pa s ml −1 for normal subjects, which match the value for the normal case at the lowest flow rate here.
Comparison of these resistance values with the idealised oral geometries of Brouns et al. (2007) shows a similar range of resistances between their cases with stenoses less than 75% and the normal case (A). The highest resistance at 30 l min −1 found here (case C, 60% constriction) is comparable with the values Brouns et al. (2007) report for cases with a degree of constriction between 85 and 90%. The differing degree of constriction causing these losses suggests that the type of idealised constriction considered by Brouns et al. (2007) does not capture the full range of loss mechanisms found in goitre cases, where the constrictions (i) tend to extend for a longer proportion of the trachea than those considered by Brouns et al. (2007), (ii) exhibit high curvature and (iii) are non axi-symmetric.
Pressure loss, which is a key measure for tracheal assessment, is not reflected purely by the minimum cross-sectional area of the trachea. Thus these results demonstrate the limitation of considering purely constrictive measures or considering stenoses to behave like orifice plates, as diagnostic tools and the need to carefully consider other geometric measures.

Mechanics of loss
The findings in the discussion above revealed similar magnitudes of pressure losses occurring in cases that had different geometric characteristics. It is necessary to characterise the mechanisms causing these losses, particularly in the cases which do not exhibit severe constriction. For this purpose, one geometry from each of the three classes (normal (case A), constricted and curved (case B) and curved, but not constricted (case E)) was analysed in  more detail to determine the loss-causing mechanisms in each of the three classes of anatomy identified. Fig. 4 shows the energy flux along the length of each trachea, at a flow rate of 30 l min −1 . The total losses in geometries B and E are approximately three times larger than in the normal case A. There is minimal loss in the normal case downstream of the first tracheal ring, whereas much of the losses in geometries B and E occur downstream of this location and are maintained from the first tracheal ring to the carina. However, there is an increase in the gradient of energy flux in geometry B in the region where the trachea expands after the constriction (around 50-60 mm), in geometry E the gradient increases at or just downstream of each peak in curvature. These increases in energy flux gradient indicate regions with higher losses and correlate to regions of high curvature and subsequent separation. Geometry C shows the highest losses throughout, particularly from the glottis to 60 mm below the first tracheal ring. Above the first ring, this is due to geometric distortion within the glottis causing multiple partial impingements of the glottal jet. The high power loss continues through the constricted section (20-35 mm) and beyond as the geometry expands and the flow separates. Only at the distal end of the trachea (>60 mm) do the losses reach an equivalent level to the other geometries.

Energy flux
The values of energy flux through the glottis and trachea shown in Fig. 4 can be put into context by comparing them to values for the work of breathing reported in the literature. Loring et al. (2009) report that in normal cases, the resistance to airflow in the airways is of order 30 mW and the total work rate is of order 300-350 mW (Milic-Emili and Petit, 1960;Loring et al., 2009). Therefore, comparison of these values with the 3 mW of energy flux required to drive flow through the normal glottis and trachea (Fig. 4, case A) reveals that approximately 10% of the airway resistance occurs in the larynx and trachea in normal cases and approximately 1% of the total work rate of breathing. The pathological cases (e.g. case C) show this tracheal contribution can increase almost 30 fold, and would therefore relate to a 30% increase in the work of breathing, all else being equal.

Flow patterns
The mechanics of loss can be explained by analysis of the flow patterns in each case. Correlation with geometric characteristics at locations of high loss allows the link between the geometric features and loss-causing mechanisms to be investigated. Therefore, the flow in three of the cases, (A, B and E) are described below. The loss-causing mechanisms in case C are similar to that in case B (although raised in magnitude), whilst case D is similar to case E.

Case A
Case A (normal geometry) exhibits a glottal jet emerging from the constriction at the glottis (a confined jet). The orientation of the glottal "throat" directs the jet towards the anterior wall where it impinges. Immediately distal to the impingement region there is breakdown of flow structure. UA MAX is elevated in this location; high velocity fluid is forced in the azimuthal direction. Further downstream the velocity profiles become blunt due to the redistribution of kinetic energy in the flow normal direction. UA MAX decays in such locations. In the distal sections of the trachea the influence of the jet is diminished and flow characteristics are related to local cross-sectional area. These observations are in accord with previous findings (Johnstone et al., 2004;Ma and Lutchen, 2006;Lin et al., 2007;Choi et al., 2009).

Case B
Case B suffers two consequences of the goitre. It is constricted as evidenced by a sustained reduction in area between x = 10 and 50 mm distal to the first tracheal ring (see Fig. 5b). Moreover, case B is displaced from the normal alignment of the trachea as shown by a gradual rise and fall in curvature between x = 20 and 65 mm (maximum curvature: 0.036 mm −1 , see Fig. 5c).
These changes alter the flow: Fig. 5a shows that the jet emanates from the glottis and impinges onto the anterior wall; the location of impingement is closely correlated with the projection of the glottal plane in its normal direction. In this upper region the influence of curvature and torsion (how sharply the centreline twists in space) seem to be negligible. Area, however, does play a role: in the constricted region where the jet impinges (10 mm < x < 50 mm), the flow accelerates as the lumen narrows. In this location, UA MEAN , shown in Fig. 5d, gradually rises as the flow stabilises and curvature begins to have an effect (see constricted region in Fig. 5b and curvature rise in Fig. 5c). Downstream of the constricted region (x > 50 mm) the lumen expands and the jet no longer dominates the flow; curvature and its direction (evidenced by the Frenet normal vectors that point towards the centre of curvature in Fig. 5g) now have a larger influence. As the flow traverses the bend, the largest force component normal to the flow is directed away from the centre of curvature (see vector superimposed onto velocity profile in Fig. 5g). This pushes the high velocity flow towards the outer wall. In this region, UA MEAN is elevated, albeit delayed due to the inertia of the flow. The reduction of UA MEAN near the end of the trachea is due to the expansion of the geometry. Finally, UA MEAN plateaus due to the moderate curvature observed in this region (x > 75 mm -dependent on flowrate). Small changes in the direction of the curvature lead to in-plane rotation of the profile, thus the multidirectional nature of curvature is important for the flow mechanics in the lower trachea.
The variation in UA MAX , shown in Fig. 5, reflects the pattern of power loss. In the jet region UA MAX is low (x < 10 mm). Following impingement and breakdown, UA MAX rises suddenly due to redistribution of high velocity fluid through the cross-section. Through the throat of the constriction UA MAX remains elevated. Anterior to the location of peak curvature, UA MAX drops as the flow detaches from the anatomical right hand side. The subsequent rise occurs due to a highly unsteady region downstream that causes power losses.
A useful indication of the state of the flow is provided by the turbulent kinetic energy (TKE), shown in Fig. 5f. TKE is initially elevated in the region of jet impingement. Progressing in the streamwise direction the TKE reduces, as is observed in normal tracheas. The main point of difference is in the expansion region (x > 50 mm) and beyond: TKE rises to a peak and then decreases gradually. The rise occurs due to flow separation (i.e. a jet like structure) from the anatomical right hand side followed by jet breakdown leading to a localised increase in TKE. TKE attenuates post this region as the energy is redistributed across the cross-section and stabilises. Although TKE essentially mimics the dynamics of UA MAX , UA MAX provides enhanced insight into the regions in which there are power losses. This is primarily because 90% of the energy flux is transported through the area defined by UA MAX .
Case E As opposed to case B, case E is only deviated from the trajectory of a normal trachea. The area, Fig. 6b, varies little along its length, while the curvature is elevated around x = 25 mm and x = 70 mm. The directions of these elevated regions are in different directions, since the goitre pushes against the trachea in the lower half of the geometry. The centre of curvature for the first peak (x = 25 mm) is on the anatomical right hand side, whilst for the second peak (x = 70 mm) it is on the anatomical left hand side.
From Fig. 6a we observe as for case B the glottal jet in the upper half of the geometry; however the alignment of the glottis means the jet does not impinge onto the wall until well downstream of the glottis (impingement occurs between planes ii and iii in Fig. 6a). As for Case B the jet dominates the flow dynamics in the upper portion of the geometry (see jet dominated location in Fig. 6a). This is captured well by UA MEAN , where the magnitude remains low for x < 45 mm and the first elevation in curvature of Fig. 6c (maximum at x = 25 mm) has little affect on the flow. Beyond this region, curvature begins to gradually influence the flow and forces the high velocity fluid from the anatomical left hand side to the right (see elevated curvature around x = 70 mm). The reasoning for this gradual influence of geometry is similar to that of case B given above. In the distal trachea, the dependence of flow on geometry is demonstrated by the elevated UA MEAN at the location of elevated curvature and the Frenet normal vectors in Fig. 6d and g, respectively.
Between x = 45 and 60 mm we observe a rapid rise in UA MAX (see Fig. 6e). The rise occurs due to redistribution of high velocity fluid in the region anterior to the impingement and breakdown of the glottal jet: post plane iii in Fig. 6a. Further downstream, UA MAX undulates due to multiple separations forming from geometric perturbations leading to localised high and low velocity regions. In the most distal part of the trachea UA MAX decays and stabilises as the geometry straightens.
TKE, shown in Fig. 6f, decays throughout the geometry, except for a small rise in the jet breakdown region. The decay prior to this region is because the jet that forms through the glottis plots a relatively stable path towards the impingement region between planes ii and iii of Fig. 6a.
The analysis of the two pathological cases highlights the importance of both UA MEAN and UA MAX : UA MEAN tracks the curvature in the lower part of the trachea; UA MAX indicates important locations of power loss. Additionally, we have highlighted that TKE fails to separate unsteady characteristics (e.g. jet oscillations) from regions with high turbulence (e.g. jet breakdown).

Significance of results
The above analysis has brought out several factors that are important in determining the overall pressure and energy losses within the trachea.
Flow in the trachea is normally dominated by the upstream laryngeal jet. Changes due to pathology influence the behaviour and breakdown of this jet in addition to changes in the affected anatomical region.
Formulae to predict tracheal losses based purely on the degree of stenosis, such as that proposed by Brouns et al. (2007) are specific to the nature of the idealisation of their geometry and thus may not represent the conditions found in vivo.
The cases presented in this study demonstrate that coupling between geometric factors such as curvature and the alignment of the laryngeal jet can be equally important as the degree of stenosis (e.g. Case E). This analysis is supported by the findings of Qi et al. (2014) who found increased pressure losses in unstenosed but deviated geometries, even with the absence of the glottis. Furthermore, in coronary stenosis assessment Escaned (2014) reports that simple stenotic measures (e.g. percentage diameter and minimal luminal area) represent an oversimplification of the geometry, thus can not characterise the problem correctly. CFD is becoming a powerful diagnostic and clinical assessment tool that has the potential to transform how patients are managed (see examples for both air flow (Hamilton et al., 2015) and blood flow (Taylor et al., 2013)). The approach presented here demonstrates that CFD can be used to determine the energy required to overcome tracheal flow resistance and its contribution to the total work of breathing based on literature values.
This study has focused on a small number of cases, due to the large computational expense of the simulations. The large cost of each simulation, combined with the need to run at different flow rates etc. mean that running a large enough volume of cases for statistical significance is not possible. The simulations in this study took an average of 3500 core hours each × three flow rates × five geometries, totalling approximately 52,500 core hours. Despite the small number of cases, they are sufficient to illustrate a gap in the current understanding of breathing losses, as the elevated losses in case E could not be predicted by consideration of constriction alone. This study has used expensive, high resolution computational simulations to identify the geometric measures and the fluidic mechanisms causing these additional losses. Now identified, future studies may be able to recreate these simulations at lower resolutions in order to simulate a larger number of cases.

Conclusions
Goitres have been used as a tool, creating suitable geometries to investigate changes in airway curvature and constriction. This study suggests that geometric factors beyond those previously considered have a profound effect on the pressure and energy losses during inspiration within the trachea. Such geometric factors include curvature, glottal alignment and deviation. In this study, the influence of pathology on these measures was found to increase the exponent of the pressure flow relationship and raise pressure and energy losses and consequently the required work of breathing.
It is convenient to divide the flow in the trachea into two regions. The first is the region in which the laryngeal jet dominates and the second is dependent on local geometry. The extent of both regions and the flow within each can be significantly influenced by pathology. As shown in this study, curvature can change the glottal jet dominated region by moving the location of impingement, thus altering the redistribution of high velocity flow. Flow in the remainder of the trachea can be influenced by curvature, constriction or by both. Curvature increases loss by forcing flow away from the centre of curvature (opposite to the direction of the Frenet normal), reducing the utilised area. Constriction is capable of inducing secondary separation regions, which further increases the loss over that due to wall frictional resistance.
Flow in all tracheas is characterised by a jet emanating from the glottis which then impinges on the tracheal wall and dissipates. However, the redistribution of the high velocity flow in pathological cases was impaired compared with normals, leading to increased power loss. Whilst the measures of utilised areas record the distribution of velocity, it has been demonstrated that they also highlight power losses and characterise internal biological flows with constrictions.
In summary it was found that (i) predictions of pressure or energy loss based purely on stenosis percentage, such as orifice plate correlations fail to capture the full range of loss mechanisms. (ii) CFD provides a tool to quantify the effort necessary to overcome tracheal flow resistance in breathing, through the reduction in energy flux. This was found to be elevated in pathological cases compared to the normal anatomy. (iii) Significant deviation enhances loss by reducing the utilised area of the flow. Quantifying this utilised area provides an intuitive tool to reduce the complex 3D flow data down to a physically meaningful 1D quantity and can also help to identify regions of enhanced power loss associated with constrictions and/or high curvature. It was found to be markedly lower in the pathological geometries compared to the normal.