Airflow through the supraglottis during inspiration

Abstract

Exercise-induced laryngeal obstruction is a paradoxical laryngeal closure during inspiration at high-intensity exercise, with supraglottic closure being most common. This study develops a model based on the computational fluid dynamics to investigate airflow velocity and pressure and the air-induced loads on the supraglottis at various inspiratory flow rates. It is found that at high flow rates, positive wall pressure is formed in the hypopharynx localise towards its lower region, while posterior supraglottic wall pressures shift from positive to negative. These findings suggest that high inspiratory flow rates may increase supraglottic pressure differentials, ultimately contributing in the collapse.

Keywords:

• Supraglottis
• inducible laryngeal obstruction
• upper respiratory tract
• inspiratory airflow
• air-induced loads

1. Introduction

Exercise-induced dyspnoea is a common ailment that impacts many patients and is often attributed to asthma or intrathoracic airway obstruction. There is an increasing amount of literature regarding collapse of the upper respiratory tract (URT) during exercise, specifically at the larynx, where a paradoxical closure can occur during inspiration (Røksund et al. 2017). This exercise-induced laryngeal obstruction (EILO) is characterised by exertional dyspnoea, inspiratory stridor, and panic reactions, which quickly dissipates after cessation of exercise (Røksund et al. 2015). Johansson et al. (2015) estimated the prevalence of EILO in the adolescent population between 5 and 7%, with the majority of those affected being female. Walsted et al. (2021) noted that this trend of female preponderance is also carried into the adult populations with diagnosed EILO.

The presentation of EILO typically occurs in the supraglottic region, but closure at the glottis, or a combination of both has been documented (Nielsen et al. 2013; Walsted et al. 2021). Supraglottic collapse during EILO involves the anteromedial movement of the arytenoid cartilages and overlying mucosa with their associated corniculate cartilages, including the posterior aspect of the aryepiglottic folds. The supraglottis is defined as the upper region of the larynx situated between the laryngeal inlet and the glottis. The mechanism behind supraglottic collapse is largely unknown but several hypotheses have been proposed. Halvorsen et al. (2017) published a statement on behalf of the European Respiratory Society and European Laryngological Society which outlines three pathophysiological mechanisms of inducible laryngeal obstruction, with mechanical insufficiency being implicated in the supraglottic obstruction observed during exercise. This hypothesis suggests that laryngeal tissue is unable to withstand forces induced by inspiratory airflow.

At present, the continuous laryngoscopy during exercise (CLE) test remains the gold standard for diagnosis and analysis of EILO. This test is not capable of studying the specific laryngeal airflow characteristics that would provide aerodynamic insights into the pathophysiology behind EILO. Walsted et al. (2018) progressed the physiological understanding of EILO through comparing the respiratory testing data of EILO patients with matched control subjects. They identified that moderate to severe symptoms were clearly observable EILO patients at peak inspiratory flows greater 180 L/min. These patients also demonstrated a higher inspiratory capacity at low exercise intensities and adopted an altered breathing pattern. Fretheim-Kelly et al. (2022) demonstrated a reliable and well tolerated method for obtaining translaryngeal pressures during the CLE test. This method offers an opportunity for additional measurements in diagnostic evaluation and treatment optimisation. However, these pressure measures do not describe pressures over the whole pressure field, or at the wall, which are valuable when attempting to deconstruct the pathophysiological mechanisms of EILO. At present, there is still a remarkable gap in knowledge about the behaviour of airflow and air-induced forces at the larynx during exercise.

Computational fluid dynamics (CFD) enables us to simulate airflow through the URT and collect data which is unattainable in the patient. Previously, CFD has been used to optimise aerosol drug delivery and predict particle deposition (Brouns et al. 2007; Shinneeb and Pollard 2012; Zhang and Kleinstreuer 2004; Koullapis et al. 2016; Cui and Gutheil 2011). These studies provide a foundation for the understanding of inspiratory fluid dynamics, but few studies have specifically investigated the laryngeal aerodynamics and air-induced forces.

Heenan et al. (2003) created three-dimensional idealised geometry of the URT based on medical imaging data and demonstrated airflow features are influenced by anatomical characteristics of the URT. This geometry identified fundamental flow features, such as flow separation, recirculation zones, and demonstrated the impact of cross-sectional area changes on airflow in the URT. Anatomical features, such as the uvula, epiglottis, and a pronounced supraglottis, were shown to influence laryngeal airflow. These findings have been supported by (Heenan et al. 2004; Ball et al. 2008a, 2008b; Pollard et al. 2012; Azarnoosh et al. 2016). Lin et al. (2007) and Choi et al. (2009) showed that complex geometry of the pharynx and larynx will impact downstream airflow, due to vortices originating in these regions and continuing towards the supraglottis and trachea. Tabe et al. (2021) have shown that the highest airflow velocities are often observed in the laryngeal jet which develops as air passes through the narrow glottis. The high-velocity laryngeal jet impinges on the anterior wall of the trachea and leads to a highly three-dimensional posterior recirculation area, as shown by Johnstone et al. (2004); Brouns et al. (2006); Collins et al. (2007); Ball et al. (2008a, 2008b); Azarnoosh et al. (2016); Pollard et al. (2012). The location of specific flow behaviours and their interaction with the URT wall is crucial for understanding to airflow-induced forces.

Respiratory flow rates vary throughout the ventilatory cycle, especially in response to changes in physical activity. Hence, it is important to consider the influence of flow rate on airflow behaviour. Koullapis et al. (2016) used a computational approach to study the air velocity profiles and particle deposition in the upper and lower respiratory tract. They demonstrated that mean flow topology was consistent across flow rates of Q = 15.2 litres per minute (L/min), 30 L/min and 60 L/min. It is known that at higher flow rates, the magnitude of the velocity profiles and turbulent kinetic energy are increased, which is supported by numerous other studies; see Kleinstreuer and Zhang (2003); Takano et al. (2006); Tabe et al. (2021); Shang et al. (2019). The changes in the wall pressures at different flow rates have been investigated by Shang et al. (2019), who found that increasing the inspiratory flow rate from 18 L/min to 50 L/min causes a significant increase in surface pressure drop at the glottis, from 60 Pa to 360 Pa, respectively. This demonstrated that pressure drop increases disproportionately to increases in flow rate. Tabe et al. (2021) experimentally identified a pressure drop of 2.5 kPa between the mouth and lower trachea at 60 L/min. However, Heenan et al. (2003) conducted laboratory experiments, which showed pressure drop, between the mouth and lower trachea, of less than 0.6 kPa at flow rates greater than 60 L/min. These dramatic differences in results are likely a consequence of geometry differences and its influence on airway resistance.

An understanding of pressure values combined with shear stresses on URT wall will provide valuable information when considering aerodynamic influences of supraglottic collapse. The influence of the flow field on these values should be considered and whether changes to URT geometry can alter the distribution of forces. Walsted et al. (2018) has shown peak inspiratory flow rates of Q = 240 L/min during exercise in EILO patients and control subjects. The majority of previous studies on inspiratory airflow utilise flow rates between 15 and 90 L/min, and do not reach values observed in maximal exercise.

In this study, a CFD model is developed to investigate the influence of several steady inspiratory flow rates on airflow and air-induced forces in the URT. We apply high inspiratory flow rates up to those observed in adults exercising at maximal intensity. This aims to develop our understanding of the aerodynamics involved in EILO, with a focus on the supraglottic and hypopharyngeal regions, as it is most prone to collapse.

In this article, we present the creation of the URT geometry and regions of interest in Section 2, followed by Section 3 on the numerical solution and its underlying theory, then the numerical set-up in Section 4. Results and discussion are given in Section 6, and this article is closed by some concluding remarks.

2. Upper respiratory tract geometry

An anonymised computed tomography scan of an adult male is used to create the three-dimensional representation of an adult URT using 3DSlicer. The geometry represents the respiratory pathway from the oral cavity to the trachea, with exclusion of the nasal cavities and nasopharynx due to the preponderance of oral breathing during maximal exercise. The inlet (mouth), outlet (trachea), and walls of the geometry are defined before exporting the geometry as a stereolithography file.

Figure 1 demonstrates the anatomical features of the URT geometry with the location and cross-sectional area of slices A-A${}^{\prime }$ to F-F${}^{\prime }$ through the supraglottic and hypopharyngeal regions. To study the wall forces acting on these regions in greater detail, they are divided into eleven regions of interest (3 hypopharyngeal and 8 supraglottic), as shown in Figure 2. The hypopharynx is defined as the region extending from the tip of the epiglottis to the laryngeal inlet, and is considered due to its upstream location to the supraglottis.

Figure 1. Upper respiratory tract geometry with key anatomical features labelled with locations and cross-sectional area of axial slices A-A${}^{\prime }$ to F-F${}^{\prime }$ indicated.

Figure 2. Location of the 11 regions of interest in the hypopharyngeal and supraglottic wall.

3. Theory and numerical solution

In this section, we discuss the theory and numerical solution of the computational fluid dynamics approach used in this study. A three dimensional, right-handed Cartesian coordinate system is used in this study, where x₁ points to the right, x₂ pointing into the paper, and x₃ points vertically up, opposite to the direction of the gravitational acceleration; see Figure 1. Indicial notation and Einstein’s summation convention are used, and subscripts after comma indicate differentiation.

The conservation of mass and momentum equations for a homogeneous, Newtonian and incompressible fluid are given by the Reynolds-Averaged Navier–Stokes (RANS) equations: (1) $${\overline{u}}_{i,i}=0,\hspace{1em}i=1,2,3$$(1) (2) $${\overline{u}}_{j,t}+{({\overline{u}}_i{\overline{u}}_j+\overline{u{\prime }_{i}u{\prime }_j})}_{,i}=g_j-\frac{1}{\rho }{\overline{p}}_{,j}+\nu {\overline{u}}_{j,ii},\hspace{1em}i,j=1,2,3$$(2) where ρ is the density of fluid, ν is the kinematic viscosity, $\overrightarrow{g}=(0,-g,0)$ is the gravitational acceleration vector, and p is the pressure. $\overrightarrow{u}=u_i\overrightarrow{e_i}$ describes the velocity field, and $\overrightarrow{e_i}$ is the unit normal vector in the i direction. $\overline{f}$ is the time-averaged value of the fluctuating arbitrary variable f.

Given the fluctuating nature of the velocities in RANS equations, we adopt the $k-\omega$ model for the closure of the equations. Given in Menter (1994), in this model, the kinematic viscosity is associated to the turbulent kinetic energy and dissipation by (3) $${\nu }_t=\frac{\omega }{\rho k},$$(3) where k is the thermal conductivity and ω is the specific turbulence dissipation rate and ν_t is the eddy-viscosity. For more details about the $k-\omega$ model used in this study; see Menter et al. (2003).

The air flow Mach number is small at all time (Ma < 0.03), and hence air is assumed incompressible in this study. The airflow is turbulent at the velocities of interest in this study, therefore it is important to consider the effect of turbulence. The k-ω model has been shown to outperform other RANS models for similar problems to that considered here, due to its treatment of the viscous effects near the wall region in wall bounded flows, similar to that in this study. Mylavarapu et al. (2009) considered an expiratory flow rate of Q = 200 L/min through the URT and demonstrated that the k-ω model showed the best agreement, against other RANS models and large eddy simulation, for static pressure measurements at the wall, when compared to experimental results. Ball et al. (2008a) and Mihaescu et al. (2008) have also found that the k-ω model shows better agreement than other RANS models to experimental results for inspiratory flows.

The RANS equations are solved numerically using a finite volume approach, which is suitable for meshes of complex geometries, such as those of the URT. Velocity-pressure coupling is achieved through the PIMPLE algorithm, a combination of the semi-implicit method for pressure-linked equations (SIMPLE) algorithm, developed by Patankar and Spalding (1972), and the pressure-implicit split-operator (PISO) algorithm, developed by Issa (1986). PIMPLE can be used to stabilise simulations of complex geometries with non-orthogonal meshes. The computations are conducted in an open source computational fluid dynamics package, namely OpenFOAM (Jasak 2009).

4. Numerical set up

Five inspiratory flow rates of Q = 30, 60, 120, 180, and 240 L/min are considered in this study. A steady volumetric flow rate is applied to the inlet boundary condition, and fixed pressure of zero at the outlet was set. In addition, the no-slip condition is applied to the URT walls. Fluid properties of ρ = 1.225 kg/m³ and ν = 1.5 × 10^{– 5}m²/s are used to represent air at 15 °C. All computational simulations are set with a variable time step to maintain a maximum Courant number of 0.5. One case is executed for 1.5 seconds of simulation, and no change in fluid flow is observed beyond 0.5 seconds. The residual convergence is set at 1 × 10^{– 5} for all variables. For the computations, two IntelXeon E5-2697A v4 processors (16 cores, 3.00 GHz) are used.

A mesh convergence study is conducted to establish the appropriate mesh for the computations. For the URT, a three-dimensional hybrid mesh consisting of hexahedral and split-hexahedral elements is generated using snappyHexMesh of OpenFOAM. Four meshes with increasing cell numbers are considered, listed in Table 1. Mesh convergence is assessed by applying an inspiratory flow rate of Q = 60 L/min and comparing the average pressure of two axial slices at the oral cavity and mid-oropharyngeal level. Table 1 shows a summary of the meshes and results. Meshes A-C were compared to Mesh D, which has an extremely high mesh density. Mesh C demonstrated the closest agreement to Mesh D with a less than 3% difference in results at the two slices and a significant reduction in computational cost; see Table 1. Therefore, we determine that Mesh C provides an accurate representation of pressure in this model and is used in the following URT simulations.

Table 1. Summary of meshes and results of convergence test. Average pressures of Slices 1 and 2 are given and the percentage difference is measured with respect to Mesh D.

Mesh ID	Cells	Slice 1 (Pa)	Differences	Slice 2 (Pa)	Differences	CPU Time (Hrs)
A	705,313	74.73	11.92 %	45.51	13.10 %	13.71
B	889,986	73.29	13.61 %	47.24	9.80 %	31.78
C	1,200,438	82.31	2.99 %	51.57	1.53 %	42.36
D	4,312,193	84.84	–	52.37	–	121.50

5. Comparisons

Results of the model given above are compared and discussed in this section. First, we present a comparison between available laboratory measurements and a computational simulation of fully-developed turbulent flow in a circular pipe to verify the numerical approach. Further, comparisons are made between the computational model of this study and published laboratory experiments and computational results of similar URT airflow cases.

5.1. Pipe flow comparison

The RANS model with k-ω turbulence closure discussed the previous sections is used to study the airflow in a pipe, and the results are compared with published laboratory measurements. Eggels et al. (1994) and Den Toonder and Nieuwstadt (1997) experimentally studied airflow through a pipe using laser Doppler velocimetry and particle image velocimetry with a Reynolds number $(Re=uD/\nu )$ of Re = 4900 and Re = 5450, respectively. In this study, a horizontal pipe with circular cross-sections with diameter D = 3 cm and length L = 20 D is considered. A uniform inlet velocity of u = 2.45 m/s is applied to produce airflow inside the pipe with Re = 4900 to match the experimental conditions.

Figure 3 shows the normalised velocity profile across the circular pipe cross-section of the experimental and computational results. The velocity magnitude is normalised by the centreline velocity u_c and the diameter is normalised by r/D, where r is the pipe radius. Outstanding agreement is observed between normalised velocity profiles across the pipe cross-section of the k-ω RANS and laboratory results of Eggels et al. (1994) and Den Toonder and Nieuwstadt (1997). The numerical approach used in this study is deemed suitable to resolve the fully developed turbulent internal flow.

Figure 3. Comparison of the results of the k-ω RANS model with laboratory experiments of Eggels et al. (1994) and Den Toonder and Nieuwstadt (1997) for velocity profile in a circular pipe.

5.2. Upper respiratory tract comparisons

To establish further confidence in the numerical model developed here, several qualitative and quantitative comparisons against experimental and computational results of flow through the URT are considered. Heenan et al. (2003) conducted an experimental study of airflow through an idealised geometry. The geometry used is modified from the one used by Stapleton et al. (2000), who constructed the URT by combining simple geometric shapes to represent the relevant anatomical features using measurements from medical imaging data; see this article for full details of the geometry and method. Heenan et al. (2003) experimentally measured pressure drop between the inlet and outlet of their URT model at a range of steady inspiratory flow rates, ranging from Q = 8 L/min to Q = 120 L/min. In our study, pressure drop values are obtained using the average pressure across the inlet and outlet for Q = 15, 30, and 60 L/min. Figure 4 plots pressure drop across the URT against flow rate. The results of the k-ω model used in our study demonstrate excellent agreement with the experimental measurements of Heenan et al. (2003).

Figure 4. Comparison of the results of the k-ω RANS model with laboratory experiments of (Heenan et al. 2003) for pressure drop between the inlet and outlet of URT at different flow rates.

Xi et al. (2018) experimentally studied the pressure drop along a realistic URT geometry, which was obtained directly from computed tomography data and provides a more detailed representation of the anatomical features. In their study, Xi et al. (2018) developed five variations of their URT geometry with different glottic areas. The glottic area of the geometry used in our study is 1.12 cm² which closely matches a geometry variation developed by Xi et al. (2018) with a glottic area of 1.14 cm² and is used for comparison. Xi et al. (2018) applied an inspiratory flow rate of Q = 15 L/min and experimentally measured pressure at nine intervals along the entire URT geometry. Our RANS k-ω model is used to generate computational comparison using the URT geometry of our study with a flow rate of Q = 15 L/min. An average gauge pressure of nine axial planes throughout the length of the URT geometry is obtained and compared to Xi et al. (2018) in Figure 5(a).

A comparison of pressure drop values between our study and Xi et al. (2018) shows near perfect agreement at the inlet and outlet, but significant variation at the larynx and the upstream regions. The variation in the two URT geometries can offer an explanation for the disparity in pressure values, but the exact cross-sectional values for the URT used in Xi et al. (2018) are unknown. The oral and pharyngeal region of our URT appears to be narrower than those in Xi et al. (2018), which leads to an increase of the fluid velocity and a subsequent decrease pressure values in this region. In Figure 5(a), both models show a sharp pressure drop at 10 cm, this pressure drop is explained by the presence of the epiglottis in the pharynx, which produces the first major constriction encountered by inspiratory airflow through the URT. At the glottis, the largest difference between pressure values of the k-ω and experimental results is observed. The cross-sectional area of the glottis between the two URT geometries are closely matched and would be expected they produce similar values. However, glottic flow field can be influenced by upstream geometries. Choi et al. (2009) and Heenan et al. (2004) demonstrate that the inter- and intra-subject variability of the URT geometry can influence the flow fields at the larynx.

Figure 5(b) shows a further comparison of the computational results of Xi et al. (2018) illustrating the similarities in the magnitude of the velocity profile along the anterior to posterior mid-point of the glottis. The RANS k-ω model show two peaks in the anterior and posterior parts of the profile, whereas the velocity profile of Xi et al. (2018) show a relatively steady increase from anterior to posterior. This is influenced by upstream geometry variation; however, as expected, the greatest values are seen posteriorly.

The comparisons of k-ω model used in this study with the laboratory experiments of Heenan et al. (2003) and Xi et al. (2018) show good overall agreement overall.

6. Results and discussion

In this section, results of the RANS $k-\omega$ model for steady airflow at rates Q = 30, 60, 120, 180, and 240 L/min through the URT are presented, with a focus on the supraglottic and hypopharyngeal regions. Results are given in two subsections; the differences in the velocity and pressure fields, and then the distribution and magnitude of the forces on the wall.

6.1. Flow fields

The velocity and pressure fields of the URT geometry are discussed here. The URT flow topology and streamlines are shown for Q = 30, 60, 120, 180, 240 L/min in Figure 6(a)–(e), respectively. The air entering the inlet passes through the oral cavity and exhibits a velocity increase, with a vortex redirecting flow from the left side of the oral cavity to the right. The average velocity of the region increases as air passes through a minor constriction to enter the oropharynx. As it continues the airflow slows as the oropharyngeal geometry expands. The airflow is directed through a constriction formed by the tip of the epiglottis, forming the epiglottic jet. This jet flow is the location of peak velocity magnitudes in the URT of 9.19, 17.89, 35.33, 53.17, 75.19 m/s were observed at inlet flow rates of Q = 30, 60, 120, 180, and 240 L/min, respectively. Airflow passing into the hypopharynx and supraglottis forms an area of recirculation anteriorly and a high velocity jet, posteriorly. In Figures 6(a)–(c), a three-dimensional vortex is present in the left of the anterior supraglottis. This vortex strengthens up to Q = 120 L/min, then appears to weaken at Q = 180 L/min and Q = 240 L/min; see Figure 6(d) and 6(e).

Figure 5. Comparison of the k-ω RANS model results with Xi et al. (2018): (a) experimental measurements for pressure drop along the URT and (b) computational results for the velocity profile at the mid-point of the glottis at Q = 15 L/min.

Figure 7 presents cross-sections of the instantaneous velocity contours and secondary streamlines in the hypopharynx and supraglottis at inlet flow rates of Q = 60 L/min and Q = 240 L/min. Labels A-F indicate the anterior aspect of the geometry, whereas labels A${}^{\prime }$-F${}^{\prime }$ indicate the posterior aspect. As the oropharyngeal flow approaches the hypopharyngeal region two counter-rotating vortices are seen in the lateral aspects, as seen in slices A-A${}^{\prime }$ of Figures 7(a) and (b). At higher flowrates, these lateral vortices in A-A${}^{\prime }$ become increasingly compact but remain in this location. Similar observations were made by Shang et al. (2019), where an increase in flow rate from Q = 15 L/min to Q = 50 L/min demonstrated variable placement of vortical structures in the laryngeal region and a weakening of supraglottic vortices. In B-B${}^{\prime }$ to F-F${}^{\prime }$ of Figures 7(a) and (b), the epiglottic jet flow persists downstream along the posterior wall of the hypopharynx and continues as the laryngeal jet into the glottic and tracheal regions. These areas of flow separation correspond to two recirculation zones, one in the supraglottis, anterior to the epiglottic jet, and another, posterior to the laryngeal jet in the trachea, downstream of the glottis; see Figure 7(b). These flow topologies are present at all flow rates studied and have been widely described in previous literature; see Liu et al. (2016); Cui and Gutheil (2011); Cui et al. (2018); Banko et al. (2015); Shang et al. (2019); Tabe et al. (2021); Zhang and Kleinstreuer (2004); Yan et al. (2019). The flow separation in these regions are due to the narrowing of the geometry leading to high-velocity flow, which when followed by a significant geometry widening with subsequent slowing of airflow, leads to complex flow patterns. The supraglottic recirculation zone shows the appearance of two counter-rotating cylindrical vortices at Q = 60 L/min but weakens and reduces to one on the left-side at Q = 240 L/min; as shown in C-C${}^{\prime }$ of Figures 7(a) and (b). In the trachea, the decay of vortices at higher flow rates have been observed in idealised (Kleinstreuer and Zhang 2003) and realistic (Shang et al. 2019; Tabe et al. 2021) geometries but seen to reduce or disappear by Q = 60 L/min. Figures 6 shows vortices in the trachea are strengthened by higher flow rates. Here we have shown that general flow topology remains the same across all flow rates considered here. However, the secondary flow features are variable as higher flow rates are reached, which can manifest as a displacement of a vortex from side to the other.

Figure 7. Mid-sagittal instantaneous velocity contours and axial cross-sections through the supraglottis with secondary velocity streamlines at (a) Q = 60 L/min and (b) Q = 240 L/min.

Shown in Figures 8(a) and (b) are the URT pressure fields for inlet flow rates of Q = 60 L/min and Q = 240 L/min, with warmer colours assigned to positive values and cooler colours representing negative values. The general distribution of pressure values are similar at different flow rates, with a significant drop in pressure is observed at the epiglottic constriction between A-A${}^{\prime }$ and B-B${}^{\prime }$ in Figures 8(a) and (b), coinciding with the epiglottic jet. Slightly downstream, significant changes in the pressure field at Q = 60 L/min and Q = 240 L/min are observed in the supraglottis and hypopharynx; see slices C-C${}^{\prime }$ to E-E${}^{\prime }$ in Figures 7(a) and (b). Slice C-C${}^{\prime }$ demonstrates the loss of positive pressure values posteriorly and a higher magnitude of negative pressure anteriorly from Q = 60 L/min to Q = 240 L/min. Slice D-D${}^{\prime }$ shows a reduction of positive pressure values in the supraglottis and localisation of positive pressure to the lower hypopharynx from Q = 60 L/min to Q = 240 L/min. Slice E-E${}^{\prime }$ shows positive pressure values in the anterolateral aspect at Q = 60 L/min, which at Q = 240 L/min are significantly reduced on the left side and lost on the right side.

6.2. Wall forces

To better understand the impact of flow rate on supraglottic collapse, the airflow-induced forces exerted on the wall are assessed. Figure 9 shows the values and distribution for gauge pressure on the surface of the URT. In general, we observe an increase in pressure values in the URT as higher flow rates are applied. The pressure distribution at different flow rates remains similar upstream of the epiglottis and in the trachea, but demonstrates variation in the supraglottic and hypopharyngeal regions. At Q = 30 L/min, positive surface pressure values continue down the posterior wall of the hypopharynx and laterally into the supraglottis, with negative values seen anteriorly coinciding with the supraglottic recirculation zone; see Figure 9(a). Whereas at higher flow rates, this region of positive surface pressure in the posterior hypopharynx becomes disconnected in the upper and middle regions and localises toward the lower hypopharynx as flow rate increases; see the middle view in Figures 9(b)–(e). This is due to the increased velocity of the epiglottic jet at higher flow rates causing a significant pressure drop in the UH and MH regions, as mentioned in Section 6.1. A consistent region of positive pressure is located at the LH region across all flow rates, while down stream of this region of the posterior wall peak negative pressure can be identified. Yan et al. (2019) computationally analysed nasal inspiration through the URT and observed the localisation of positive pressures in the lower hypopharyngeal region at an inspiratory flow rate of Q = 90 L/min, which was not present at lower flow rates Q = 50 L/min and Q = 30 L/min.

Figure 8. Mid-sagittal pressure field contours and axial cross-sections through the supraglottis with secondary streamlines at (a) Q = 60 L/min and (b) Q = 240 L/min.

Figure 9. Pressure distribution on the URT wall for each flowrate presented from three views with left, middle, and right images being the left lateral, posterior, and right lateral views, respectively. (a) 30 L/min, (b) 60 L/min, (c) 120 L/min, (d) 180 L/min, and (e) 240 L/min.

In the context of EILO, this demonstrates regions of the high positive and low negative pressures, creating a pressure differential, in close proximity to the regions prone to collapse which is reinforced as flow rate is increased. In the lateral supraglottis, the increase in flow rate reduces the pressure, which continues inferiorly to leave small areas of positive pressure near the glottis. Downstream of the glottis, an area of positive pressure on the anterior wall is observed due to the high-velocity flow of the laryngeal jet interacting with the anterior trachea. Posterior to this, in the region of the tracheal recirculation zone, negative gauge pressures is observed.

Figure 10(a) shows the pressure force acting on the regions of interest presented in Figure 2. The pressure force is obtained by $F=\int p\mathrm{d}A$, where p is the pressure, and A is the surface area. The pressure forces in all hypopharyngeal subdivisions are positive at Q = 30 L/min then become negative at higher flow rates Q > 120 L/min in the UH and MH. The LH region exhibits a positive normal force at all flow rates in the LH, which peaks at 180 L/min with an average pressure of 84 Pa over the region. The high-pressure observed in the LH, occurs due to collision of the epiglottic jet against the downstream wall before it is redirected anteriorly through laryngeal inlet. The maximal negative forces were observed in the ULA and URA regions at all flow rates, shown in Figure 10(a). These regions include the lateral margins of the epiglottic tip where peak velocities are observed in Figure 6. In Figure 10(a), the URP and ULP regions demonstrate inconsistency in their normal force values as higher flow rates are applied to the inlet. The URP has relatively low positive normal force at flow rates up to Q = 180 L/min, but at Q = 240 L/min drops significantly and becomes negative. The normal force of the ULP region is positive at low flow rates, which become negative at Q > 120 L/min and reaches a peak negative value at Q = 180 L/min.

Figure 10. Comparison of regions of interest. (a) Pressure forces acting on regions of interest in the original URT geometry and (b) normalised force comparison on regions of interest with different epiglottic angles. The forces are normalised with respect to the values of the lower hypopharyngeal (LH) region at each configuration. LH, lower hypopharynx; MH, middle hypopharynx; UH, upper hypopharynx; ULA, upper left anterior; URA, upper right anterior; ULP, upper left posterior; URP, upper right posterior; LLA, lower left anterior; LRA, lower right anterior; LLP. lower left posterior; lower right posterior; see Figure 2 for the relative positions.

Figure 6. Lateral (left) and posterior (right) views of URT streamlines coloured by velocity magnitude for each inspiratory flow rate. (a) 30 L/min, (b) 60 L/min, (c) 120 L/min, (d) 180 L/min, and (e) 240 L/min.

Wall shear stress (WSS) was also considered at all flow rates. The region with the highest WSS values across all regions at all flow rates was UH with a peak value of 1.78 × 10^{– 3} N at 240 L/min, much smaller than the pressure force, and is not remarkable. The contribution of shear stresses on the loads exerted by air flow on the wall is very small in comparison to the normal forces and are unlikely to contribute to the collapse of these regions.

6.3. Effect of geometry

To investigate the significance of the geometry on the results discussed in this study, we consider two modified geometries in this section. The two geometric modifications from the original are achieved through anterior displacement of the epiglottis and are characterised by the angle of the epiglottis relative to the posterior oropharyngeal wall. The original geometry exhibits an epiglottic angle of 32 degrees, whereas modification 1 and modification 2 have angles of 27 degrees and 22 degrees, respectively. As previously discussed, the epiglottic projection into the URT causes a substantial increase in flow velocity, which likely dictates the downstream flow features observed in the supraglottis. Therefore, using these geometric modifications, we investigate effect of a geometric change on the force profile over the supraglottic and hypopharyngeal regions of interest. These cases were run with a steady, uniform flow rate of 180 L/min at the inlet. As previously described in Section 6.2, the LH region demonstrates the largest positive force across all regions and flow rates. Therefore, to better understand the distribution of forces and establish a profile of wall forces, each region is normalised by presenting it as a ratio of the LH force values for the respective geometry.

Figure 10(b) shows the distribution of forces relative to the LH values with close agreement between the 27- and 22-degree modifications. In the 32-degree geometry, the regions near the epiglottis exhibit much lower forces. The regions of the lower supraglottis (LLA, LRA, LLP, and LRP) across the modifications exhibit close agreement in values which suggests that the magnitudes of these forces are consistent, relative to the LH values of the same geometric configuration in this model. Relative forces in the upper regions of the supraglottis and hypopharynx are more variable, where the regions closest to the epiglottis (UH, MH, ULA, and URA) demonstrate the greatest variation. This is unremarkable as these are the regions in close proximity to the geometric modifications. This is in contrast with the posterior areas, URP and ULP, which are more closely related. These results demonstrate that wall forces in the lower supraglottic regions relative to the LH region are consistent between geometries of varying epiglottic geometries.

6.4. Further discussion

Our results demonstrate agreement with the general flow features identified by previous airflow studies through a URT geometry (Heenan et al. 2003, 2004; Ball et al. 2008a,2008b; Pollard et al. 2012; Azarnoosh et al. 2016). It is shown that the flow topology remains consistent across all flow rates, but secondary flow features exhibit variation in their location as flow rates increases within the same geometry. The pressure field and subsequent wall forces exhibit variation in magnitude and distribution in the posterior supraglottic and lower hypopharyngeal regions, which is clearly evident at flow rates of 180 L/min and greater. In the context of EILO, these regions are in close proximity to areas that are observed to collapse during high intensity exercise. The presence of the epiglottic jet appears to cause high positive pressure forces in the lower hypopharyngeal region that increase with flow rate. This coupled with the negative pressure forces seen at higher flow rates in the upper posterior regions of the laryngeal inlet may be suggestive of a force differential facilitating the aerodynamic collapse of the supraglottis. This is in-line with the findings of Walsted et al. (2018) who identified that all EILO patients in their cohort with an inspiratory flow rate of 180 L/min, or higher, at mid-inspiratory volume had moderate to severe laryngeal obstruction.

The flow rates considered in this computational model were steady and uniformly applied at the inlet.

We also observed that the velocity and pressure fields exhibit maximal values in close proximity to the epiglottis. The modifications considered here showed that the relative force profiles demonstrate variations in regions related to the epiglottis, but regions downstream were relatively consistent between the geometric configurations. This suggests that changes to the geometry may influence the force distribution locally but do not drastically influence relative wall forces and load distribution further away from the site of geometric change.

In this study, a single URT geometry (with some variable features) with a steady inspiratory flow is considered. The geometry omits several specific anatomical features that were not retained during creation of the geometry, which led to a simplification of the anatomy when compared to that observed in the patient. Future work should incorporate more detailed laryngeal anatomy and various anatomical configurations into the geometry to observe their influence on airflow. In real life, the breathing cycle demonstrates acceleration and deceleration of flow velocity occurring in inspiratory and expiratory directions. The unsteady oscillating flow pattern may differ in characteristics from the steady unidirectional flow studied in this work and should be explored in future research to ascertain its potential effect on the aerodynamic supraglottic collapse.

7. Conclusions

A k-ω RANS model is developed to numerically investigate the impact of steady inspiratory flow rates on airflow and air-induced forces in a computational URT model. The model is verified in comparison with available laboratory experiments and other numerical data. The URT geometry is developed from computed tomography data and several inspiratory flow rates are considered, ranging from Q = 30 L/min to Q = 240 L/min, to represent values observed during maximal exercise.

In this study, we find that the general flow topology and distribution of surface forces remain largely similar across the URT at lower flow rates, 30 L/min and 60 L/min. However, the supraglottic region exhibits variation in the velocity and pressure fields at flow rates greater than 120 L/min, which consequently leads to altered distribution of wall forces. These wall forces remain positive and localise in the lower hypopharynx, while becoming increasingly negative at the posterior supraglottis. These findings suggest that a hypopharyngeal-supraglottic pressure differential, exacerbated by high flow rates, may be a contributing factor that leads to the supraglottic collapse seen in EILO patients. This study provides the first analysis of inspiratory airflow at rates greater than 180 L/min. Further investigation is required to study the influence of geometric configurations and time-varying flow seen during the breathing cycle.

Disclosure statement

The authors report there are no competing interests to declare.

Funding

The author(s) reported there is no funding associated with the work featured in this article.