Inflow boundary profile prescription for numerical simulation of nasal airflow

D. J. Taylor, D. J. Doorly, R. C. Schroter

Abstract

Knowledge of how air flows through the nasal passages relies heavily on model studies, as the complexity and relative inaccessibility of the anatomy prevents detailed in vivo measurement. Almost all models to date fail to incorporate the geometry of the external nose, instead employing a truncated inflow. Typically, flow is specified to enter the model domain either directly at the nares (nostrils), or via an artificial pipe inflow tract attached to the nares. This study investigates the effect of the inflow geometry on flow predictions during steady nasal inspiration. Models that fully replicate the internal and external nasal airways of two anatomically distinct subjects are used as a reference to compare the effects of common inflow treatments on physiologically relevant quantities including regional wall shear stress and particle residence time distributions. Inflow geometry truncation is found to affect flow predictions significantly, though slightly less so for the subject displaying more pronounced passage area contraction up to the internal nasal valve. For both subject geometries, a tapered pipe inflow provides a better approximation to the natural inflow than a blunt velocity profile applied to the nares. Computational modelling issues are also briefly outlined, by comparing quantities predicted using different surface tessellations, and by evaluation of domain-splitting techniques.

1. Introduction

Computational modelling of physiological flows in realistic anatomical conduits has developed rapidly in the last two decades. Developments in the fields of imaging, computational methods and computer technology now enable flow predictions to be obtained in patient-specific geometries. There has been an enormous growth in the number of computational studies addressing cardiovascular flows in particular, and, to a growing extent, also respiratory flows. In both application areas, the complexity of the flow system is too great to be fully modelled at present. A commonly applied means to reduce computational effort is to exclude the geometry outside the region of interest, with truncation applied at the inflow and outflow boundaries of the modelled domain.

For investigations of nasal inspiration, most studies (with a few exceptions, such as Croce et al. 2006) ignore the external nose, assuming that during inspiration the inhaled air is imposed either directly at the nares (nostrils), or that it is supplied directly to the nares via an artificial inflow tract. It is clearly preferable to incorporate an adequate definition of the external nose when modelling nasal inspiration. Although improved computational power now makes this possible, nevertheless, the effects of artificially truncated inflows need to be considered for nasal inspiration (as indeed for other physiological flows) for two reasons. Firstly, estimates of the modelling error associated with truncation (derived from comparison with full modelling) provide a guide to determine whether the extra computational cost is needed to attain the level of fidelity required. Secondly, and arguably a more important need, is the fact that many in vitro model studies have been, and are still, restricted to truncated inflows. Thus, although errors can be avoided in computational modelling, reliance is still placed on model studies for estimates of quantities such as aerosol deposition that are difficult to predict. Comparative studies are thus needed to reconcile newer computational models of the full domain, with experiments in truncated models.

This paper aims to address the effect of inflow boundary profile prescription on the patterns of inspired air by simulating laminar airflow through two healthy unilateral nasal cavities. The following sections briefly detail the physiology of the nose and the background to previous model studies, before describing the methods used and the results obtained. While clearly focused on modelling of the nasal airways, the approach employed and the methods used may be of interest for application to other physiological flow studies.

2. Background

The human nose performs three physiological functions, namely: (i) heating and humidification; (ii) olfaction; and (iii) filtration and defence. Though each task separately constitutes a difficult challenge, the nose attains a remarkable level of performance in all three. Firstly, inspired air largely attains body temperature and saturation before reaching the glottis (Mygind & Dahl 1998; Wolf et al. 2004). Secondly, although olfactory sensitivity depends on a host of factors, including the particular odorant molecule, threshold values for olfactory detection in the low tens of parts per billion have been reported (Walker et al. 2003). Finally, all particles larger than 5 µm are removed and approximately 50 per cent of particles between 2 and 4 µm are captured by the nasal mucosa (Wolf et al. 2004).

In engineering terms, much of the nasal function is achieved by passive means, whereby the control of the airflow to accomplish the physiological tasks is largely effected by the architecture or form of the internal passages. However, there is inherent conflict in the modes of airflow required to achieve these tasks: to promote efficient heating and humidification calls for rapid flow transit with minimal pressure loss, whereas to facilitate sampling and mass transfer of odorant molecules to the enervated epithelium calls for more lengthy transit times. Moreover, filtration with minimal pressure loss is inherently difficult. The rather convoluted form of the nasal airways is a response to these diverse challenges and is discussed elsewhere (e.g. Lang 1989; Wolf et al. 2004). However, although the overall principles are believed to be known, the mechanisms by which the airflow and transport processes operate are as yet poorly understood.

Understanding the mechanisms of airflow is not only relevant to physiology, but may influence fields such as surgical planning, drug delivery and toxicology. In vivo investigation of nasal airflow is limited by the inaccessibility and complexity of the nasal passages, reducing assessment to gross measures of nasal function. Techniques such as rhinomanometry and acoustic rhinometry (Hilberg 2002) can be used to identify gross nasal obstruction but are incapable of providing details of the nature of the obstruction or its impact on nasal performance.

The nature of in vivo airflow can be analysed in far more detail using both in vitro and in silico simulation techniques. Procedures have now been developed to produce high-fidelity transparent models, which replicate the nasal anatomy through rapid prototyping and inverse casting procedures (e.g. Hopkins et al. 2000; Kelly et al. 2000; Kim & Chung 2004; Taylor et al. 2005; Doorly et al. 2008a). While physical models can provide detailed measures of deposition, pressure drop and velocity distributions, they involve a large commitment of resources. Computational simulations offer a potentially less demanding option, but care needs to be taken to ensure that the computational measures can be relied upon.

Previously, computational simulations have been applied by, among others, Keyhani et al. (1997); Naftali et al. (1998); Zhao et al. (2004); Elad et al. (2008) to investigate processes such as olfaction and air-conditioning. Other studies using computational simulations aim to address a documented need for a more objective assessment of nasal function (Eccles 2000; Hanif et al. 2000) and to provide surgeons with a pre-operative insight into the post-operative implications of different surgical procedures (Lindemann et al. 2005; Wexler et al. 2005) as well as consideration of pathological conditions (Garcia et al. 2007). Computational simulations can also provide probabilistic maps of particle deposition in the nasal passages (Kleven et al. 2005; Schroeter et al. 2006; Shi et al. 2006, 2007). Such information could also benefit toxicologists by identifying how inhaled airborne pathogens or particulates are transported to the nasal passageways from the external environment.

The computational results presented here are from a study in which anatomically accurate, twice scale silicone phantoms were fabricated, of the right nasal airways of two healthy subjects. By controlling and balancing the model scale, fluid viscosity and flow velocity, the experiment can be made to replicate the dynamics of the real flow, though with the benefits of higher spatial resolution (using a larger scale) and temporal resolution (using liquid instead of air effectively slows the motion by an order of magnitude). The resultant flow-field will be kinematically identical to that of a 1:1 scale model with circulating air, as discussed in Doorly et al. (2008a); a comprehensive account of scaling is given in Barenblatt (2003).

The fabrication process comprised the initial translation of in vivo images of the airways into physical rapid-prototypes, and inverse casting to yield the final silicone replicas, as detailed in Taylor et al. (2005), Taylor (2009) and Doorly et al. (2008a). To ensure that computational simulations could be optimally compared with future experimental studies, the computational wall meshes were generated via high-resolution medical imaging of these silicone phantoms. Such an approach has the advantage of allowing any observed differences between the airflow behaviour in computational and experimental studies to be characterized, free from geometric ambiguity associated with the errors in the translation of the computational anatomy to a model replica. Thus, discrepancies can be attributed to computational modelling inaccuracies, which may arise due to, for example, poorly prescribed boundary conditions or insufficient mesh resolution.

3. Methods

Two computed tomography (CT) image datasets, derived by retrospective examination of a large set of patient CT image records, were obtained with permission from the ENT Department of St Mary's Hospital, London. These were segmented (Amira™, Mercury Computer Systems Inc.) by manual delineation of the nasal boundaries on successive images to yield a three-dimensional surface triangulation for each image dataset. Their respective image acquisition parameters are shown in table 1. Openings (ostia) to the sinuses were omitted as they presented minimal cross-sectional area and were deemed to have a negligible impact on the gross airflow patterns. In-house bi-Laplacian smoothing techniques were implemented to remove pixellation artefacts and to provide a more realistic, smooth surface finish to the geometries, as discussed by Gambaruto et al. (2008). Inflow and outflow extensions were attached to the replica geometry definitions at this stage to allow manufactured physical models to be connected to a flow loop for experimental studies.

View this table:
Table 1.

Properties of the computed tomography (CT) image acquisition of two nominally normal nasal anatomies. The pertinent imaging parameters of the original CT acquisition are shown, while the parameters of CT acquisition of the silicone phantoms (which are double scale) are shown in parentheses.

Twice scale unilateral silicone phantoms were created of both subject anatomies using rapid-prototype and inverse-casting procedures as described elsewhere (Taylor et al. 2005; Doorly et al. 2008a). The twice scale silicone phantoms were scanned using high-resolution CT, with image parameters as given in table 1. In-house segmentation procedures were used to automatically reconstruct a computer definition of the surface boundaries of each silicone model. With no extraneous ostia, and given the high signal-to-noise levels in contrast differences between the silicone- and air-filled portions of the model, automatic methods proved insensitive to the segmentation threshold level. The surfaces of the original and scanned silicone phantom were registered using iterative closest point procedures and it was found that the root mean square distance between their surfaces were less than 0.29 mm at 1:1 scale (i.e. less than an imaging pixel).

The surface meshes were manually sectioned normal to an expected nominal mean flow direction using a CAD program (Rhinoceros, Robert McNeel & Associates), yielding the area distribution curves (similar to acoustic rhinometry measurements) and sections of the anatomy shown in figure 1a,c, respectively. The symbols dotted along these curves denote the locations of slices extracted through both geometries, which are shown in figure 1c. The nasal airway of subject B is revealed as having a much larger volume, although the characteristic forms of the two subject area distributions show good correlation, as illustrated by alignments of the peaks and troughs of cross-sectional area distributions of both subject anatomies.

Figure 1.

Anatomical description of the right nasal passageways for subjects A and B. (a) Cross-sectional area (CSA) distributions, extracted normal to an expected mean flow path of both nasal geometries are shown; X represents the inter-slice distance between successive area-centroids. The symbols denote the locations of the slices shown in (c), and are shown as surface lines in figure 8. Subject A, dashed line with filled square; subject B, solid line with filled circle. (b) Magnified views of two slices extracted through the volume meshes used for in silico simulations depict typical mesh resolution at the apex of the nasal cavity and on a slice through the anterior tip of the middle turbinate for subjects A and B (top and bottom, respectively); extracted slices (3 and 6) are depicted in (c) and denoted using vertical arrows in figure 8. (c) Selected slices extracted through both geometries show the calibre of the passageway of subject A as being much smaller than that of subject B; the slice orientations are shown in figure 8 and are consistently ordered 1–9 from nostril to nasopharynx. The slicing locations are depicted with symbols on the area distributions of (a).

These CT datasets are broadly representative of each end of the normal healthy range, as indicated in table 2. Specifically, the surface area to volume ratio (SAVR, cf. Garcia et al. 2007), defined as the ratio of nasal airway surface area to volume for the region extending from the nostril until the end of the septum, was 0.92 and 0.53 mm−1, respectively, for subjects A and B. The difference in calibre between the subjects is reflected by the fact that subject A has a passage volume much smaller than that of subject B, yet they have similar surface areas. Hence, the SAVR of subject A is nearly twice that of subject B. Subjects A and B also display differing nasal valve orientations and nasal vestibule area distributions, which is significant given that the nasal valve typically accounts for 50 per cent of the resistance of the entire respiratory tract (Wolf et al. 2004). Moreover, gross flow patterns within the cavity depend strongly on the angular orientation and rate of contraction of flow up to the nasal valve.

View this table:
Table 2.

Nasal passage dimensions for both subject anatomies. Typical ranges of nostril and nasal valve cross-sectional areas are 50–130 × 10−6 and 20–60 × 10−6 m2, respectively (Lang 1989). The Reynolds numbers are based on the hydraulic diameter (4 × area/perimeter) of the nasal valve and for constant inspiration at 100 ml s−1.

4. Mesh generation

The surface meshes, reconstructed from the scanned silicone phantoms, were used for computational fluid dynamics simulations to provide an optimal geometry match with future experimental studies. Commercial CAD software (Rhinoceros) was used to produce the surface meshes for each simulation, as shown in figure 2, namely: (a) configurations (i) and (ii) consisting of the nasal passageway from the nostril to nasopharynx; (b) configurations (iii) and (iv), which include a convergent pipe inflow tract; (c) configuration (v), which includes an external face, also segmented from the CT datasets, and bounded by a box extending approximately three face widths, two face depths and one nose length from the face. As the neck is typically omitted from CT image acquisitions, due to the susceptibility of the thyroid gland to ionizing radiation, the lower boundary of the reconstructed face (or chin) was extended downward to mimic the neck.

Figure 2.

Different inflow boundary profile configurations are illustrated using subject A, though the same configurations were applied to subject B. These are representative of commonly encountered in vivo, in vitro and in silico configurations, where: (a) (i) a flat (constant) velocity profile was applied normal to the nostril; (ii) the prescribed velocity profile at the nostril was extracted from a separate, prior, simulation of flow through a convergent pipe inflow up to the nasal valve; (b) (iii) and (iv) depict the application of a flat and parabolic velocity profile, respectively, at the entrance to a convergent pipe inflow; and (c) (v) depicts the application of pressure boundary conditions to a simulation including the external face, allowing flow to develop naturally towards the nostril (shown using trajectories of passive particles).

It was found that the flow-field about the external face acts as a sink-like flow, with air drawn uniformly towards the nostril. The magnitude of the velocity field decays significantly with distance from the nostril (by an order of magnitude within 1 cm of the nostril, as discussed in Doorly et al. 2008b). Hence, physiologically realistic inflow conditions can be generated by simulation of a limited portion of the external flow-field (thus reducing the computational cost). Although not incorporated in this research, computational modelling could be simplified by the introduction of a symmetry boundary condition (bisecting the face sagittally) to further reduce the computational effort required to resolve the flow-field about the external face.

The triangular surface tessellations were re-meshed (Gambit v. 2.3.16, Fluent Inc., Lebanon, USA) to provide high-quality surface elements for computational simulations. Volume meshing (TGrid v. 4.0.16, Fluent Inc., Lebanon, USA) of each configuration yielded an unstructured mesh with a core of tetrahedral elements and five prism layers adjacent to the wall, for improved boundary layer resolution.

The solutions were performed on medium- to high- resolution meshes (approx. 12–24 million prism and tetrahedron elements), with the first prism height typically 1 per cent of local passageway gap width; the typical densities of volume meshes are indicated in figure 1b. The average wall shear stress (WSS) and overall pressure drop for subject A were found to be under-predicted by approximately 0.9 and 0.5 per cent, respectively, when comparing a mesh density of 15 million cells (which is typical of those presented in this article) to one of 29 million cells, as discussed in more detail by Taylor (2009), so that the simulations may be considered nominally mesh-independent at these levels of grid density. Furthermore, it was found that a simulation comprising approximately 3.4 million elements under-predicted the same measures by approximately 3 and 2.9 per cent, respectively, which may be adequate for many applications. Simulations conducted on a computational mesh comprising approximately 3.4 million elements typically require approximately 7 h to converge (below 10−7 for scaled residuals) when run in parallel on two 2.66 GHz quad core processors and require approximately 4 Gb of RAM.

5. Computational fluid dynamics simulations

The discretized Navier–Stokes equations, which govern fluid flow, were solved on the unstructured meshes using a commercially available finite volume solver (Fluent v. 6.2.16, Fluent Inc., Lebanon, USA). Assuming incompressible flow of a Newtonian fluid and ignoring body force terms, the momentum equation reduces toEmbedded Image 5.1where v, ρ and μ represent, respectively, the fluid velocity, density and dynamic viscosity applied to the momentum equation. Neglecting temperature effects and density variations, the continuity equation,Embedded Image 5.2provides the zero divergence constraint.

Simulations of steady flow were run in double precision using a three-dimensional laminar viscous solver with the third-order MUSCL (Monotone Upstream-Centred Schemes for Conservation Laws) interpolation scheme to improve the spatial accuracy in solving the momentum equation. The SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithm was used to enforce mass conservation and to obtain the pressure field using second-order interpolation. The continuity and momentum equations were solved sequentially using a segregated implicit solver to reduce memory requirements and preserve convergence stability, while node-based gradients were used to calculate derivatives for the unstructured meshes.

The nasal walls were modelled as rigid, no-slip boundaries. The Newtonian fluid used was air with a viscosity and density set to 1.7894 × 10−5 kg m−1 s−1 and 1.225 kg m−3, respectively. All simulations were conducted at a constant volumetric flow rate of 100 ml s−1, representative of quiet restful inspiration. For simulations inclusive of a bounded external face, pressure inlet and outlet boundary conditions were constrained to target the desired mass flow rate. In all other simulations, the mass flow rate was enforced by the application of a velocity profile at the inlet (Dirichlet) in combination with a pressure outlet boundary condition.

For each subject geometry computations were undertaken of configurations representative of those typically encountered in studies of nasal airflow (cf. figure 2) listed as follows.

(i) A flat (alternatively referred to as constant or plug) velocity profile was applied at the nostril, as used in several previous studies due to its simplicity.

(ii) A quasi-developed profile was applied at the nostril, which was extracted from a separate prior simulation of flow in a pipe inflow up to the internal nasal valve. Inspiratory flow through the nasal valve is normally unidirectional, and the parabolic character of attached internal flows implies that the extent of the upstream influence is small. Although splitting the domain into two overlapped regions increases the total computational cost, the memory requirements for either of the split computations are reduced.

(iii) A blunt and (iv) a parabolic profile, respectively, were applied at the entrance of a pipe inflow tract. This is of interest to experimentalists as it addresses concerns regarding the required entrance length for developed flow, and the effect of the flow loop profile (at the entrance to the convergent pipe inflow) on experimentally resolved airflow patterns.

(v) Pressure boundary conditions were applied at the inlet and exit of the flow domain to allow the flow to develop around the external face and enter the nostril at its naturally developed angle, representing nasal inspiration in its physiological configuration. For realism, an equal flux was allowed to develop at each nostril, while the contralateral nasal valve presented a second outflow, the internal flow in the corresponding nasal passageway was not computed.

6. Results

The predicted WSS distributions, shown in figure 3, are a useful measure in assessing intra-subject sensitivity to inflow condition. Not only is WSS a sensitive indicator of mechanical processes, such as heat and mass exchange, but as Elad et al. (2006) point out, may exert mechano-transductional effects on the nasal epithelium. Contralateral views depicting flooded contour plots of WSS magnitude for subjects A and B are shown in figure 3a,b, respectively. The differences in WSS distribution for simulations including the pipe inflow tract have been omitted as they are similar to those of the flat and face simulations.

Figure 3.

Contours of the magnitude of WSS (WSS − τ*), normalized by the mean WSS magnitude (table 4) as viewed from each side, are shown for (a) subject A and (b) subject B. The simulations for a flat velocity profile at the nostril (top of a and b) and including the external face (bottom of a and b) correlate strongly, though noticeable differences occur in the region of the olfactory cleft and inferior meatus. Surface streamlines of WSS (i.e. streamlines following the WSS vector along the nasal boundary) indicate the underlying flow direction and patterns; note the stagnation of the inspiratory jet in the roof of the cavity, the large anterior recirculation region forward of this and above the jet issuing from the nasal valve. Transverse WSS streamlines adjacent to the nasal valve of subject B depict the detachment of the jet issuing from the nasal valve from the lateral wall of the nasal passageway. Broken lines denote the locations of slices (3 and 6) referred to in figures 1 and 4.

The magnitude of the WSS is highly sensitive to the internal flow structure and reveals subtle changes of flow patterns according to the applied boundary profile. The WSS is also of physiological interest as it is a strong indicator of the regional exchange burden of the nasal mucosa. Thus, the high levels of WSS on the cavity walls, both in the vicinity of the internal nasal valve and on either side of the anterior head of the middle turbinate (which is hidden in the views shown), indicate thin boundary layers that facilitate efficient heat and mass transfer; in contrast, the reduced shear stresses in the olfactory cleft indicate longer residence times for odorant retention and benefits transport to the olfactory epithelium.

For subject A, the surface shear lines shown in figure 3a indicate flow features such as the location of the stagnation of the inspiratory jet in the roof of the cavity and the large anterior flow recirculation region. The application of a flat velocity profile normal to the nostril plane is seen to lead to increased frictional losses indicated by the increased levels of WSS in the nasal vestibule. The global surface streamline and WSS patterns show considerable similarity, irrespective of the prescribed inflow configuration (with pockets of high and low WSS localized to consistent regions). However, some local differences are apparent, principally in regions of reduced flow, such as in the apex of the cavity and inferior meatus.

For subject B, the relatively increased passageway calibre compared to subject A results in reduced mean airflow speeds and frictional losses. The mean WSS magnitude for subject B (cf. table 4) is approximately three times lower than that for subject A, although the overall patterns of WSS are broadly similar in both the subjects. Again, the application of a flat velocity profile at the tip of the naris induces high WSS levels in the nasal vestibule of subject B, and the reattachment location shows notable anterior movement when compared with the configuration inclusive of the external face. The upper and lower, left-most views of the WSS distributions of subject B in figure 3b reveal detachment of the inspiratory jet from the right lateral wall of the nasal cavity, just downstream of the nasal valve. The surface streamlines can be seen to converge, indicating departure of flow from the wall. In addition, the WSS in the inferior meatus of subject B is greatly reduced due to its high patency and low flow rate as the core flow, carried by the jet issuing from the nasal valve, bypasses its entrance. Though the recirculating velocities in this region are weak, there are marked differences in inferior meatal recirculation patterns for the two boundary profiles—notably, the posterior reattachment in this region forms an expansive region of slow retrograde flow with application of a flat inflow profile.

The changes in internal flow pattern for the different inflow configurations are also highlighted by examination of the WSS distributions along the boundary of two slices extracted through the computational wall boundary (indicated as slices 3 and 6 in figure 3 and arrowed in figure 8, which is discussed later). The first slice is extracted through the anterior section of the nasal cavity, approximately midway between the nasal valve and anterior head of the middle turbinate, and the second in the mid-length of the cavity, intersecting the region of the olfactory cleft. The WSS profiles extracted at these slice locations are plotted in figure 4, where the horizontal axis (P*) represents the normalized perimeter. These correlate very strongly for subject A, but there are evident differences for subject B, corresponding to variations in inflow; this is consistent with the mean normalized differences shown in table 3.

Figure 4.

Plots of normalized WSS magnitude (τ*) along the perimeter of the slices extracted through the computational wall meshes are shown for (a) subject A and (b) subject B. The slice locations correspond to those identified in figure 1 (slices 3 and 6) and are denoted in figure 8 by vertical arrows. Distance (P*) has been normalized by perimeter length and WSS (τ*) by μŪ/dh at the nasal valve. The form of each extracted slice is shown to the right, where coloured symbols are used to represent the unwound perimeter location also shown beneath the horizontal axis of the plots to the left. Again, for subject A, little excursion is observed for differing inflow profiles, though subject B continues to show strong changes. Blue line, flat; red line, pipe inflow; black line, face.

View this table:
Table 3.

The mean (Embedded Image) and s.d. (σΘ*) of the absolute differences between the WSS magnitudes for different inflow configurations are shown for each subject (expressed as a percentage of the local mean values). The mean deviations for the three configurations incorporating a pipe inflow show negligible change, ranging from 2 to 6% for subject A and 6 to 11% for subject B.

View this table:
Table 4.

An assessment of the global measures (pressure drop (ΔP), mean residence times (Embedded Image) and average WSS (Embedded Image) across the nasal passageway) shows limited variation with inflow profile; however, regional measures (integrated Embedded Image and average (Embedded Image) WSS and volumetric flux (Q) to the region of the olfactory cleft) highlight significant regional augmentation of flow pattern, particularly for subject A. The volume flux into the upper regions of the nasal cavity (Q) measured as a percentage of the volumetric flow rate through the nose (100 × 10−6 m3 s−1) are tabulated for the varying inflow profiles. Hence, for both subjects 10 per cent of inspired air reached the region containing the olfactory cleft for a flat inflow velocity profile at the nostril. The mean residence times (Embedded Image) for 40 000 passive particles tracked through the computational domain—from nostril to outflow—are also tabulated. (Because little variation was evident between pipe inflow configurations—(ii), (iii) and (iv)—these are tabulated as an average of the three.)

Similarities between the pipe inflow and face configurations for both the subjects are shown in figure 5 by means of flooded contours of velocity magnitude extracted at the nasal valve, together with streamlines indicating in-plane motion. The imposition of a flat velocity profile at the nostril forces high-speed air through the narrow apex of the nasal valve and results in a differing overall distribution of velocity. The pipe inflow and face configurations show broad similarity, with peak flow biased more towards the inferior portion of the nasal valve. This provides an indication that the resultant flow patterns within the main nasal cavity will differ more significantly with the application of a blunt velocity profile than for the pipe inflow configuration, when compared with the most realistic configuration that includes the external face.

Figure 5.

Profiles of velocity magnitude normalized by mean velocity (U*) are shown for subjects A (left) and B (right), and for each inflow configuration. The streamlines indicate in-plane velocity while the flooded contours represent the normalized magnitude of velocity. Again, the calibre of subject B can be seen to be much greater than that of subject A. For both subjects the application of a flat (vertical) velocity profile at the nostril forces fluid into the apex of the nasal valve, in marked contrast to the pipe inflow and face configurations, where peak flow is located to the posterior and inferior portions of the nostril and nasal valve, respectively.

The pressure drop through each anatomy, plotted as pressure drop coefficient (Cp—the local pressure normalized by ½ρŪ2 at the nasal valve), is shown in figure 6. However, note that the absolute pressure loss of subject A is approximately four times that of subject B (see summary of table 4, which is discussed later). The normalized pressure was calculated using the average of the area-weighted total pressure across sequential planar slices extracted through the computational volume (shown in figure 1 and as surface lines in figure 8). The horizontal axis represents the distance along the area-centroids of the extracted slices. This overall pressure drop represents the contribution of the nasal airways to the respiratory workload and can be monitored clinically (e.g. by rhinomanometry) as a global means of predicting nasal function. Prescription of a flat inflow velocity profile is shown to exaggerate the initial pressure drop, as expected, since unrealistic WSSs within the nasal vestibule are induced by the initial boundary layer. However, particularly significant features of these results are that: (i) the pressure drop curves for each subject show negligible difference downstream of the nasal valve and (ii) the pressure drop for subject A is nearly the same for face and pipe inflow configurations.

Figure 6.

Pressure drop coefficient (Cp; left) and particle residence time distributions (right) are shown for (a) subject A and (b) subject B. Cp was calculated using the average area-weighted total pressure across slices extracted through each volume mesh, and for different inflow boundary profiles (circles, flat; squares, pipe inflow; triangles, face; the slice locations are depicted in figures 1 and 8). The distances between successive area-centroids (X) are plotted against pressure coefficient, which is normalized by the dynamic pressure at the nasal valve (Cp = ΔPρŪ2). Posterior to the nasal valve the pressure curves for all configurations drop in unison, though deviations are more noticeable for subject B; the pressure drop curves, through the passageway of subject A, are nearly identical when comparing the pipe inflow and face configurations. Probability mass functions (PMFs) are used to depict the residence time (t*, normalized by mean residence times shown in table 4) distributions of 40 000 passive particles tracked through the computational domain. The convergence of residence time distributions is depicted for subject A using a flat profile; the inset plot represents the profile for clouds of 40 000 and 20 000 particles (left and right bars, respectively). Blue bars, flat; red bars, pipe inflow; black bars, face.

The residence time distributions, also shown in figure 6, are depicted as probability mass functions (PMFs), which were calculated using the transit times of 40 000 passive particles; these were initially uniformly distributed at the nostril and tracked until exiting the nasal passageway. The overall mean transit times show some differences according to how the inflow is specified (quantitative measures are summarized in table 4). Intra-subject differences can be related to changes of flow pattern within the cavity; for example depending on the partitioning of flow more particles may be entrained in regions of slow or retrograde flow. The residence time distributions for subject A are tightly grouped, for all inflow profiles. However, the forms of these distributions are not as consistent for subject B, particularly for low transit times, indicating greater sensitivity of the convective transport by the flow to the inflow condition; this effect was preserved irrespective of the selected transit time ranges used to compute the distributions. For example, in the case of subject B, approximately three to four times the number of ‘fast’ particles transit the domain within approximately 0.12 s (corresponding to a normalized time t* = 0.5 or half the mean transit time) for a pipe inflow compared with a flat or face inflow, while little change in transit rates for any of the four transit time ranges of particles is observed for subject A.

While changing patterns of WSS corresponding to different inflow specifications can be observed visually (cf. figure 3), procedures to extract quantitative measures of variation are desirable given the importance of WSS. Since different inflow configurations resulted in different surface triangulations, even for the same individual, direct differencing of nodal quantities cannot be applied. Here, to derive the mean and variance of differences in WSS patterns, inverse distance interpolation (appendix A) was used to map the WSS magnitudes between the wall meshes of each subject, thus allowing for the fact that the different inflow configurations resulted in dissimilar surface triangulations even for the same individual. (Such an approach may prove useful in other areas of physiological flow modelling.)

Relative variations of the WSS magnitude (Θ*) were then calculated as a percentage of the absolute difference between the WSS magnitudes of each mesh compared, normalized by their mean local value (appendix A). This allows changes to be weighted by their local significance rather than an absolute global measure. As shown in figure 7, comparisons for each subject highlight in particular: (i) notable differences in the region of the olfactory cleft, in both cases, where typical deviations were found to be in excess of 50 per cent; (ii) greater deviations between flat and face inflows (upper comparison for each subject) than between pipe and face inflows (lower comparison for each subject); and (iii) greater overall deviations for subject B than for subject A. Noteworthy also are deviations in the inferior meatus, caused by gross changes of flow pattern, which can be seen in the WSS streamlines of figure 3.

Figure 7.

Colour contours of the relative deviation of the WSS magnitude (Θ*), expressed as a percentage of the mean local value, on the outer lateral walls of the cavity are shown for (a) subject A and (b) subject B. For each subject the WSS values obtained by the application of a flat profile at the nostril (top) and pipe inflow (bottom) are subtracted from the WSS values from the simulation including the external face and normalized using their respective mean values (see appendix A for details). The region of the olfactory cleft shows notable deviations, as does the inferior meatus. When comparing the flat and face profiles (top of a and b), the spatial extent of the deviations is considerably larger than is found when comparing the pipe inflow and face profiles.

The area-weighted means and s.d.'s of the intra-individual differences are tabulated in table 3 (see appendix A for further details on the mathematical implementation). The deviations are tightly grouped for subject A, with mean deviations between 24 and 28 per cent reflecting the fact that all profiles show similar differences. For subject B, deviations are considerably larger ranging between 31 and 44 per cent. Little discernible difference was found between the mean deviations for the three configurations incorporating a pipe inflow (i.e. simulations (ii), (iii) and (iv) shown in figure 2); these ranged between 2 and 6 per cent for subject A and 6 and 11 per cent for subject B. For both anatomies the pipe inflow correlates most strongly with the face configuration.

The observed deviations of WSS can be attributed to changes in internal flow patterns, and were found to be markedly different in the roof of the nasal cavity, as shown in figure 7. Though the olfactory cleft is subjacent to the cribiform plate, bounded laterally by the superior turbinate and nasal septum and limited coronally by the sphenoidal sinus (Lang 1989), the distribution of enervated olfactory neuroepithelium is not restricted to the olfactory cleft; its spatial extent varies widely in the literature, e.g. from 1–3 cm2 (Leopold et al. 2000) to as high as 10 cm2 (Illum 2003; Bear et al. 2006). The current study identifies the ‘olfactory region’ as the region shaded orange in figure 8, which includes the dorsal aspects of the roof of the cavity (both the olfactory cleft and the region directly surrounding it; this accounts for approximately 9 cm2 of the nasal surface area for both the subjects).

Figure 8.

The two nasal anatomies are shown at left for (a) subject A and (b) subject B with the region containing the olfactory cleft coloured orange. Slices extracted through the computational domain are shown traced around the outer wall of each anatomy and are referred to in figures 1 and 6. Probability mass functions (PMFs) are used to depict the relative extent of exposure of the olfactory mucosa to differing ranges of WSS magnitude (τ); these are tabulated in table 4. Blue bars, flat; red bars, pipe inflow; black bars, face.

The volume flux (Q) to this region and the associated regional forces on the nasal mucosa (Embedded Image and Embedded Image) were calculated for each boundary profile prescription and are tabulated in table 4. Even relatively minor alterations of the trajectory of the inspiratory jet, and hence its stagnation on the middle turbinate, can dramatically bias the volume flux penetrating into this region. It was found, for subject A, that the volume flux to the olfactory region was lowest for a flat profile. This was opposite in the case of subject B which showed an increased flux for this configuration (table 4). It should be noted that volume flux to this region does not necessarily correlate directly to the regional distribution of WSS on the nasal mucosa; if the inspiratory jet entering the olfactory region is spatially diffuse, the integrated or average WSS can increase even for a decreasing volume flux (table 4). This is particularly noticeable for subject B, with a 20 per cent decrease in flux to the olfactory region resulting in an increased average and integrated WSS across the olfactory epithelium.

Variations of flow with inflow profile are categorized in table 4 as either global measures (pressure drop, mean transit time and an average of the overall WSS) or measures localized to the olfactory region (integrated and average WSS, and volume flux). These global measures show little variability with changes of inflow profile, particularly for subject A, although regional changes may be important. For example, the average and integrated WSS in the olfactory region doubles for subject A (with a 50% increase in flux to the olfactory region) when comparing the flat to the pipe inflow configurations. Subject B shows less dependence on boundary prescription for regional measures but exhibits stronger global variations.

PMFs, shown in figure 8, were calculated for both subject anatomies; these show the distribution of WSS across the wall of the olfactory region (coloured orange). For subject A the columns are closely grouped except for low WSS for the flat profile. Subject B shows larger overall variation of WSS distribution with applied boundary profile. Interestingly, the application of a flat profile exposes more of the olfactory mucosa to low WSS than for either a pipe inflow or face configuration. This is true of both subjects irrespective of the olfactory fluxes which are higher for pipe inflow and face configurations of subject A and lower for subject B. For instance, approximately three times the olfactory wall area is exposed to the highest levels of WSS (30 × 10−3 Pa and above) for subject A for the pipe inflow compared with that of the flat velocity applied at the nostril. This is partially balanced for the lowest levels of WSS (6 × 10−3 Pa and below), which have approximately one-third of the area exposed for the pipe inflow compared with that of the flat profile.

7. Discussion

The computational simulations presented here represent a range of inflow configurations typically encountered in model studies of airflow in the human nasal cavity, namely the application of: (a) a flat velocity profile at the nostril; (b) a flat and parabolic profile at the entrance to a convergent pipe inflow to the nose; and (c) pressure boundary conditions to a simulation including the external face. The influence of the inflow configuration used on the resolved flow patterns is investigated within the right nasal passageways of two healthy subjects. The study is intended to assess the potential scale of error in various approximations to the physiological state (i.e. the natural nasal inflow), both with reference to the many existing studies, which have employed artificial inflow conditions, and for guidance in future studies in which a natural inflow cannot be incorporated into the modelling methods.

To reduce the resources required for computational simulations it has been shown (using configuration (ii), cf. figure 2) that the nasal valve provides a convenient site for splitting an otherwise oversized mesh into manageable blocks. This would also allow a velocity profile, extracted from a limited prior simulation, to be applied to multiple simulations of flow through similar anatomical geometries, e.g. for simulation of different surgical procedures within the same anatomy (turbinate thinning, sinus surgery, correction of septal deviation, etc.).

Computational simulations using configurations (iii) and (iv) highlight the negligible impact of applying either a blunt or parabolic boundary profile at the entrance of the convergent inflow tract. This is helpful to experimentalists as it negates the need for flow loops extensive enough to provide a fully developed (or known) velocity profile at the replica model inflow (which may be considerable for scaled-up models).

For configuration (v), which incorporates the external face, inspired air is allowed to develop naturally up to the nostril. Apart from being physiologically the most realistic configuration, it is of potential interest to toxicologists allowing the identification of the spatial origins of inhaled particulates or noxious gases. It is also noteworthy that the airflow patterns correlated more strongly when comparing the pipe inflow configuration (rather than flat profile) to that of the face configuration, for both subject anatomies.

It has been shown that for gross flow measures, such as pressure drop and residence time distributions, the inflow profile has negligible effect for subject A (which has the narrowest calibre and provides the strongest contraction through the nasal vestibule) but exhibits notable variation for subject B. For instance, the pressure drop for subject A is nearly identical for pipe inflow and face configurations; even for a flat profile, if the non-physiological pressure drop through the nasal vestibule is accounted for. Croce et al. (2006) in their model study found that entrance effects had little influence on the pressure drop measured through their model, noting that it was largely predetermined by the degree of contraction of the nasal vestibule. With configuration (i), i.e. the application of a flat velocity profile at the nostril, subjects A and B exhibit similar relative pressure losses of 24 and 18 per cent up to the nasal valve and 64 and 60 per cent up to the anterior head of the middle turbinate, respectively. This is within the range of pressure loss measured by Croce et al. (2006) and Schreck et al. (1993) for an equivalent section of the airway up to the head of the anterior turbinate who, respectively, quote a pressure loss of 76 and 43 per cent.

The distributions of WSS were used to quantify the influence of the inflow velocity profile on the characteristics of airflow within the nasal passageways. Methods were introduced which account for the localized impact of airflow variations with differing inflow configuration; subtle changes of global airflow patterns can result in significant effects in regions of low flow (such as the olfactory cleft where local airflow characteristics influence olfaction). WSS patterns were found to exhibit local deviations of over 50 per cent and highlighted that deviations (for both subjects) were clearly stronger when comparing a flat inflow profile (rather than a pipe inflow) to a simulation including the external face. The most notable deviations were found in the region of the olfactory cleft.

Given the importance of the local flow environment to olfaction, probability mass functions of WSS were used to further illustrate the dependence of the boundary configuration on flow in the region of the olfactory cleft. The most significant intra-subject variations were found when comparing differing inflow profiles for subject A, while subject B showed moderately less dependence. For instance, the volume flux to the region of the olfactory cleft, for subject A, increased by 50 per cent when comparing a flat profile with that of the pipe inflow. For a flat velocity profile 10 per cent of inspired air was found to reach the region of the olfactory cleft for both the subjects.

These findings were found to be in broad agreement with the computational simulations conducted by Keyhani et al. (1995) who found that typically 10 per cent of inspired air reached the olfactory cleft, and was nominally independent of flow rate (125–200 ml s−1). Likewise, Hahn et al. (1993) in their experimental study of the same anatomy found that approximately 14 per cent of the nasal flow was conducted through the olfactory cleft, independently of the flow rate (180–1100 ml s−1). Similarly, Zhao et al. (2004) determined olfactory flux to be 2 and 8 per cent, respectively, for the left and right nasal passageways, and found that this was strongly dependent on the anatomy of the olfactory cleft.

Previous studies in multiple nasal anatomies and for flow rates up to 200 ml s−1 (in a single cavity) found that airflow within the main nasal passageway was predominantly laminar, though increasing instability at the margins of the inspiratory jet was observed (Doorly et al. 2008b). The patterns of airflow within the main cavity were markedly similar across a range of flow rates (100–200 ml s−1), which suggest that a comparable level of sensitivity (of measures such as WSS) to inflow configuration would be maintained for these flow rates. However, when extrapolating to the highest flow rates consideration must be given to the potential for displacement of the alar wall (the side wall of the nostril), as this could significantly modify the geometry of the nasal vestibule.

8. Conclusions

The use of a truncated inflow in models of nasal inspiration affects the flow within the nasal cavity, despite the area contraction up to the internal nasal valve. The level of sensitivity to inflow prescription was found to be relatively low for qualitative flow patterns and gross flow measures. Specifically, variations of the order of 10 per cent for pressure drop, integrated WSS and even mean residence times were observed; the degree of variation was found to be reduced for the case of a more constricted nasal vestibule. However, high levels of sensitivity (of up to 100%) to inflow conditions were observed for measures restricted to the region of the olfactory cleft, such as regional WSS and olfactory flux. Hence, where the regional flow patterns are of importance (e.g. for prediction of odorant transport to the olfactory epithelium) careful consideration must be given to the prescribed inflow conditions.

Acknowledgements

This research was funded by the Biotechnology and Biological Sciences Research Council (E18557 and BB/E02344/1). The authors also wish to extend their appreciation to Dr K. L. Lee, Dr A. Gambaruto, Dr P. Franke, Dr V. Franke for their help and advice which contributed to this research. We also appreciate the assistance of the ENT Department, in particular Dr R. Almeyda and Mr N. Tolley, and the Department of Radiology, St Mary's Hospital, London.

Appendix A. Methods

As the surface triangulations for different inflow configurations are not the same, inverse distance interpolation was used to directly compare the WSS distributions for different boundary configurations. The resolution of the computational meshes was found to be sufficiently fine to allow accurate interpolation applying inverse distance weighting using the Shepard method (Shepard 1968) as shown in equation (A 1):Embedded Image A 1where PT and Pi are the locations of the mesh node on the target and source meshes, respectively. In this case interpolation was performed between each node on the target mesh to the nearest element on the source mesh (N = 3) and distance was weighted linearly (X = 1).

Relative variations of the WSS magnitude (Θ*) were then calculated using equation (A 2) as a percentage of the absolute difference between the WSS magnitudes of the target mesh (τT) and each corresponding interpolated value (τI), normalized by their mean local value:Embedded Image A 2where τT and τI represent the WSS magnitude at each mesh node of the target mesh and their respective interpolated values from each element on the source mesh, respectively.

Area-weighting was also incorporated to remove the bias which would have resulted from non-uniformly sized mesh triangulations. Equations (A 3) and (A 4) were used to calculate the mean (Embedded Image) and standard deviation (σΘ*) of the percentage differences of WSS magnitude calculated using equation (A 2):Embedded Image A 3Embedded Image A 4where ai represents the area of each vertex (one-third of the cumulative area of adjacent elements) and Θ*i represents the relative deviation of WSS magnitude expressed as a percentage of the local mean for each of the N wall vertices.

Footnotes

    • Received July 21, 2009.
    • Accepted August 7, 2009.

References

View Abstract