Abstract

Current approaches to model nasal airflow are reviewed in this study, and new findings presented. These new results make use of improvements to computational and experimental techniques and resources, which now allow key dynamical features to be investigated, and offer rational procedures to relate variations in anatomical form. Specifically, both replica and simplified airways of a single subject were investigated and compared with the replica airways of two other individuals with overtly differing geometries. Procedures to characterize and compare complex nasal airway geometry are first outlined. It is then shown that coupled computational and experimental studies, capable of obtaining highly resolved data, reveal internal flow structures in both intrinsically steady and unsteady situations. The results presented demonstrate that the intimate relation between nasal form and flow can be explored in greater detail than hitherto possible. By outlining means to compare complex airway geometries and demonstrating the effects of rational geometric simplification on the flow structure, this work offers a fresh approach to studies of how natural conduits guide and control flow. The concepts and tools address issues that are thus generic to flow studies in other physiological systems.

Keywords:

1. Introduction

The relation between anatomical form and function is of enduring interest, being central to our understanding of physiological systems, and providing lessons for engineering design. The interaction between an organism and its environment involves processes such as heat and mass transfer, and the mechanisms employed to effect these processes are constrained to be feasible by the governing physical laws. The impact of these constraints on biological form is strikingly evident when respiratory function is considered. Over short distances, diffusion alone is an effective mechanism for oxygen exchange, but the transport time rises with distance squared. The quantitative basis for this ‘diffusion limitation’ was described in Krogh's (1919) original work. Consequently, species with characteristic body radius beyond approximately 1 mm rely on convective mechanics to transport quantities within a diffusion-effective proximity.

The effectiveness of convective transport depends on the flow dynamics, in turn characterized by the Reynolds number (Re) relating the inertial/viscous force balance, i.e.where ρ is the fluid density; μ is the dynamic viscosity; Q is the flow rate; U is the mean velocity; and d is the hydraulic diameter (4×cross-sectional area/wetted perimeter) of the conduit. Where Re<1 (e.g. in the smallest airways), momentum transport by diffusion rapidly permeates the flow, so the viscous stresses ensure that the flow pattern is tightly controlled by the local conduit boundaries. Conversely, where Re≫1, inertial dominance lessens local boundary control of the flow. Complex patterns such as detached vortices or zones of flow reversal may occur, and correspondingly uneven distributions of viscous stresses; at sufficiently high Re, spontaneous flow irregularity and eventually turbulence arise.

Flow in the nasal airways can attain peak values of Re in excess of a thousand and is thus inertially dominated. Nasal flow is further complicated by the intricate form of the cavity and further compounded by the diversity of subject anatomy. Studies of nasal airflow thus illustrate many of the difficulties posed by complex flow and form alike; while in this work attention is focused on a single conduit in the body, the techniques and approach to modelling described are relevant to other physiological flows.

The paper is organized as follows. The relevant anatomy is first described and flow modelling reviewed with an example to illustrate. The representation of anatomical form is then considered in more detail. In the succeeding sections, two key flow features, namely the impact of the inspiratory jet on the middle turbinate and the limits of steady laminar flow, are discussed. Rather than confining attention to a single nasal geometry, comparisons are made with two distinct replica geometries. High-speed particle image velocimetry (PIV) illustrates flow beyond the limits of truly steady behaviour, where direct computation is still challenging both to undertake reliably and to interpret. Finally, new perspectives on modelling and scope for future directions are summarized.

2. Nasal airway form and the modelling of function

The anatomy, physiology and function of the nose are described in several works, including Proctor (1982, 1986) and Mygind & Dahl (1998). As part of the respiratory system, the nasal airways warm, humidify and cleanse inspired air, to protect the delicate tissues of the lungs. Measurements and estimated overall performance figures in these roles are given in Wolf et al. (2004). Defence mechanisms, which include alterations in passage size, and mucosal secretions, highlight the adaptability of the nasal airways to respond to external challenges. As an olfactory organ the nose samples inspired air, directing a portion to sensory receptors where it is temporarily retained to facilitate odorant molecule capture. That the nose accomplishes such disparate tasks is remarkable, since efficient air conditioning requires a rapid throughput, whereas chemical sensing benefits from longer time scales. The nasal airways also occupy limited space, measuring approximately 10 cm front to back and with a main cavity height of approximately 5 cm top to bottom.

Unsurprisingly, the passageways of the internal nose are intricate, narrow and difficult to access, effectively precluding detailed in vivo measurement using current technology. The overall pressure loss sustained by flow through the nasal cavity can be determined in vivo, but is a poor determinant of healthy nasal function, furnishing no regional information. Sensitive measures cannot be deduced from the correlation of overall quantities with anatomical form alone: given Re>O(100) within the cavity, flow is not necessarily evenly distributed in the airspace. Other means are needed to obtain detailed information on the flow, and modelling presently affords the only viable approach. Given the marked inter-species variation in airway morphology, human-specific modelling is essential and little reliance can be placed on the translation of animal model findings.

Nasal airway modelling is of use not only for respiratory physiology but also for current and putative applications. These include planning and assessment of surgical interventions (e.g. Garcia et al. 2007) and the design of precision-metered drug delivery systems to target the systemic blood circulation via the highly vascularized nasal mucosa.

The first task in modelling nasal airflow is to define the bounding geometry. Various models have previously been used, some based on anatomical idealizations (Naftali et al. 2005) and others on semi-idealizations (Hörschler et al. 2003) or anatomical replicas. Replica geometries have been created ex vivo, by casting or plastination of cadavers (e.g. Croce et al. 2006) or derived from in vivo computed tomography (CT) or magnetic resonance imaging (MRI) data (e.g. Hahn et al. 1993; Hopkins et al. 2000; Hausserman et al. 2002; Kim & Chung 2004). Each approach has respective advantages and disadvantages; while cadaveric geometries can be imaged or cast to provide well-defined boundaries, their airways may display artificially high degrees of decongestion, as the mucosal tissue shrinks during post-mortem specimen preparation.

Although in principle in vivo imaging records actual geometry, the resolution achievable is limited, particularly with MRI. The resulting uncertainty in airway boundary data has hitherto received scant attention in flow modelling. A complicating feature is that nasal geometry is never fully fixed; for example, natural cyclic variation in airway calibre associated with alternating states of congestion/decongestion is reported to occur in a significant proportion of the population (Eccles 2000); furthermore, responses to irritants from the external environment can cause changes. Nonetheless, the uncertainty in replicating a specific anatomy should be distinguished from the real variations in geometry.

This work relies on pre-existing patient CT data, obtained with permission from the ear, nose and throat (ENT) surgical department at St Mary's Hospital, Paddington, London. In a small proportion of patients referred for CT scanning of airways or sinuses, the cavity anatomy is subsequently deemed normal by consultant radiologist or surgeon. Such data provide a basis for modelling, though is generally of inferior resolution vis-à-vis cadaveric data. For example, in vivo CT data are typically obtained via axial slices (i.e. in the transverse or foot-to-head direction) to limit radiation exposure, whereas a coronal acquisition (i.e. in the anterior-to-posterior direction) would provide better definition of fine passage structures. The original dataset comprises 82 axially acquired images captured to a 512×512 pixel matrix. Overlapped image slices were acquired, with slice thickness of 1.3 mm and a spacing of 0.7 mm, resulting in an in-plane pixel size of 0.39×0.39 mm. Newer multi-slice CT scanners could now allow better resolution: 0.5 mm slice thickness and 0.2×0.2 mm in-plane, thus reducing boundary uncertainty considerably.

Figure 1a shows a slice through a coronal re-projection of CT data; figure 1b shows a computer reconstruction of a part of the right nasal cavity superimposed on a sagittal (i.e. orthogonal to axial and coronal planes) slice through the CT data. The distribution of computed wall shear stress (WSS) on the right nasal cavity wall is indicated by colour-coded contours, as will be discussed later; the left nasal cavity and septal wall are removed to allow an unobstructed view of the inferior and middle turbinates. The slice location is marked with a vertical white bar in figure 1b, and as shown lies roughly mid-way through the cavity.

The main anatomical features are identified in figure 1a,b. The septum (s) divides the cavity into left and right sides of approximately equal volumes, from the nares (twin anterior openings at the nostrils) to the nasal choanae. These twin posterior openings funnel inspired air to the nasopharynx (np). In the cavity on each side of the septum, the passageway is much reduced by the turbinates (labelled mt and it), curled bony structures covered in erectile tissue that project into the airway. Below each turbinate runs the correspondingly named passageway or meatus, e.g. the inferior and middle meatuses. The superior turbinate is smaller and not readily discernible in figure 1a, being located towards the rear (posterior) of the cavity. Also evident in the CT image are the air-filled sinuses; the left maxillary sinus is labelled (ms), and a complex network of ethmoidal air spaces is visible on each side of the upper cavity.

Figure 1c shows the variation in passage area along the inspiratory flow path from its entrance at the nares to the nasopharynx, for two different individuals: the curve labelled R-1 corresponds to the subject 1, as illustrated in figure 1a,b, while the curve labelled R-3 is for a different subject. Passageway cross sections for the two individuals at three, approximately corresponding locations are shown in figure 1d, where the left member of each pair corresponds to geometry R-1. The CT slice of figure 1a is approximately mid-way between the locations (2) and (3) marked in figure 1c.

On both sides of the nose, between the nasal vestibule and the main cavity, the airway narrows and is normally the smallest at the internal nasal valve (nv), marked by the inclined white bar in figure 1b, and arrows labelled 1 in figure 1c. For subject R-1, the internal nasal valve measures 42 mm2 in area at life scale. The reported area range in normals is variously quoted as 20–60 mm2 (Lang 1989), 30–40 mm2 (Proctor 1982) and 46–115 mm2 by Cakmak et al. (2003); the last of these used more recent CT image data.

The flow computational results incorporated in figure 1b reveal much of the relation between the flow and the dominant geometric features, and are worth briefly considering before describing the modelling techniques in more detail. The computation is for inspiratory flow at a steady rate of 100 ml s−1 (i.e. a minimal level of quiet, restful breathing). The nasal valve of subject R-1 displays a pronounced narrowing to a small opening, resulting in a pronounced cavity jet inflow. The imprint of this jet on the cavity walls is reflected in the shear stress levels (values in Pa), indicated by surface colouring. High values occur near the minimum area at the nasal valve and where flow strikes the middle turbinate. Low shear stress reveals low flow in the upper reaches of the cavity, with the superimposed black surface streamlines showing recirculation; flow entering the main cavity thus does not instantly expand to fill the airway uniformly. WSS is not merely a sensitive indicator of flow, and key determinant of heat and mass transfer, but provides a mechanotransduction mechanism for the epithelial cells lining the cavity to respond to the external environment, as discussed by Elad et al. (2006). Their reported maximum values of WSS correspond well with those observed here (0.5 Pa), bearing in mind the somewhat higher flow rate they considered.

The white arrows also superimposed in figure 1b indicate the area-weighted mean flow path, defined as the locus of the flux-weighted centroid of successive sections arranged transverse to the flow. Determining the centroid and orientation of the sections requires limited iteration, with the section orientation continually adjusted until it is orthogonal to the local average flux. In figure 1c, the variation in total cross-sectional passage area along a mean flow path of two different anatomies is compared, specifically models R-1 and R-3. The minimum on each area curve corresponds to the internal nasal valve; at 96 mm2, the area for model R-3 is more than twice that of R-1, also the degree of flow contraction from naris to internal valve is proportionately less for model R-3.

The disparity in passage area between subjects reflects different degrees of congestion, with R-3 highly decongested. For model R-1, the total pressure drop from nares to nasopharynx is 9.7 Pa at a flow rate of 100 ml s−1. In this case, 24% of the total pressure drop is accumulated by the flow in reaching the internal nasal valve. Pressure then falls rapidly as flow transits the nasal valve, with 64% of the total occurring by the anterior head of the middle turbinate. Thereafter, pressure drops more slowly, with 76% lost by the mid-meatal section indicated. By contrast, for R-3 at the same inspiratory rate, the total pressure drop is only 2.0 Pa, though the distribution of loss is similar: 18% occurring by the point of minimum cross-sectional area, 60% by the anterior head of the middle turbinate and 75% by a corresponding mid-meatal section.

Returning to the modelling procedures, translating images to a computational representation of the airway is the first task. Despite the enormous density contrast between air and tissue, the airway boundary could not be delineated automatically, and semi-manual image segmentation, using commercial software (Amira, Mercury Computer Systems, Inc., UK), was applied. Firstly, it was necessary to exclude openings (ostia) to the sinuses to restrict the model domain. Secondly, fine dividing structures, particularly in the upper cavity, suffer from partial volume effects; by reviewing the data in axial and coronal projections, an operator familiar with the anatomy can check consistency and manually exclude false connections apparent at constant threshold levels. Little flow penetrates to the upper cavity, so inaccurate boundary definition in the region of the ethmoidal cells and superior meatus is unlikely to alter gross flow characteristics. Indeed for model R-1, the short superior meatus was deliberately excluded. By contrast, the apparent patency of the lower meatus at this location depends on image threshold level (figure 1a(ii)). For these data, operators with different anatomical knowledge (one, a trainee ENT surgeon) produced slightly different segmentations: in one case a fully patent and in the other a partly patent lower meatus. Higher resolution CT would reduce image uncertainty, though some degree of operator variability may remain an issue.

Segmentation defines a pixellated (stair-stepping) airway boundary that must be smoothed to provide a realistic definition of the bounding surface. In-house smoothing algorithms, developed using the methods of Kobbelt et al. (1998) and incorporating curvature-adapted smoothing, were applied as described by Gambaruto et al. (2008).

Once the airway bounding surface is defined, computational- or experimental-based modelling can be applied. Computational modelling used the commercial package TGrid (Ansys) to generate a hybrid mesh, comprising a five-layered prismatic mesh adjacent to the walls, and a graded tetrahedral interior mesh. The height of the first prismatic layer was typically 0.04 mm (within 2% of the local passage width), with a minimum of 0.6 million elements for the surface and 15.5 million for the interior; a tiny portion is shown in figure 1a(iii). Mesh convergence checks employed up to 29 million interior elements to fill a single half-cavity. The element count is far higher than hitherto used, but ensures very high spatial resolution of both WSS and its gradients; other convergence studies in steady laminar flow (Franke et al. (2005), and to be reported elsewhere) show that lower resolutions (3–5 million elements in a carefully graded, hybrid mesh) are adequate to resolve pressure, velocity and overall shear stress.

3. Variation in nasal airway form and the modelling of complex geometry

Techniques similar to those mentioned above, transforming image data to a computational model and prediction of flow, have been used by Zhao et al. (2004), Schroeter et al. (2006) and Shi et al. (2006) to investigate the transport and deposition of particles, such as fine aerosols, and odorant gases. The issue of variability in airway geometry has however scarcely been addressed: current imaging admits some uncertainty in airway definition, and, more significantly, there are large natural inter- and intra-individual variations. Consequently, variations between the replica models and the effect of model idealization should be investigated. Comprehensive population-based studies would require prodigious effort for image segmentation and flow solution; the aim here is to describe suitable methodologies and identify how characteristic geometric features affect the flow.

To examine sensitivity to perturbations of a single geometry, three alternative models for subject 1 were constructed: model R-1a by a trainee ENT surgeon, familiar with the anatomy, model R-1b, by a less experienced researcher, and the third, model R-1c, comprising a slight simplification of R-1a using Fourier descriptors by a second researcher (figure 2a). Of these nominally identical replicas, most studies used R-1a, referred to elsewhere in the text simply as R-1. The major differences between R-1a and R-1b are: firstly, R-1b displayed a partly blocked lower meatus (instead of fully patent); secondly, the nasopharynx of R-1b expanded to full width, i.e. it included that of the opposite left nasal airway, instead of being split along an assumed plane of symmetry. A further, idealized representation of model 1, generated using ellipses or combination of ellipses to simplify the cross section, is shown labelled as I-1. Note that for the idealized model, the septum has been straightened and the nasal valve simplified, but the cross-sectional area was adjusted to replicate that of R-1, taking care to preserve the calibre of the olfactory cleft.

Variation in nasal cavity geometry. (a) Models labelled R-1a, R-1b and R-1c corresponding to slightly different replica airways of subject 1 derived by different users. Model I-1 is an idealization using elliptical fits to contour boundary, while local cross-sectional area is preserved. (b) Mid-section of replica R-1a removed to show airway ‘skeleton’. Using skeletonization, topological form and local airway calibre of models R-1a, R-2 and R-3 are combined to produce averaged airway (arrows). Inset shows mid-cavity sections; note averaging reduces local variations, see straightened septum. (c) Lower meatus (lm) skeleton of model R-3 and its curvature distribution κ(s). POD analysis extracts dominant modes from large dataset; reconstruction using dominant modes shows that these may be used to characterize form.

In order to compare different nasal anatomies, compact representations of complex geometries are needed; one approach is based on separate extraction of topological form, the skeleton, and the local passage width or calibre as described next. Skeletonization in two dimensions has been introduced (Blum 1967) to extract and describe shape, and since then used largely in image processing and machine vision (Jain 1988; Sonka 1995; Weeks 1996). Methods to describe the skeleton of an object include the crossing points of concentric fronts expanding from the object border (Blum), the locus of centres of locally maximally inscribed discs or spheres and level sets based on distance maps (Kimmel et al. 1995).

For this study, skeletonization was applied to the stack of two-dimensional coronal slices defining the boundary contour (Gambaruto 2007), and the essential steps are as follows.

For each coronal slice, the centroid C of the area enclosed by the bounding airway contour Γ is found; Γ is defined by point pairs, e.g. {(xi, zi)}.

Orthogonal axes are defined, with one axis aligned in the direction of the point on Γ furthest from C, say (xf, zf), together with an associated uniform sampling grid. The grid spacing is typically set to allow 1000 increments along the major axis.

At each point on Γ, the outward normal ni is computed, starting with the furthest point (i=f), where the direction is unambiguous and then marching around Γ.

The vector dP from each sampling grid point P to the nearest point Q on Γ is found, and the signed distance, d(P)=dP.nq, calculated from the scalar product of dP and the normal, where the sign discriminates between interior and exterior points, and d(P)=0 if P is on the boundary. (Note d(P) is analogous to an elevation map, with the curve along the ridge joining peak values corresponding to the skeleton that is to be extracted.)

A simple 3×3 mask filter applied to the map d(P) yields a cluster of points about the skeleton, which are thinned and connected to form curves. The procedure also allows an automatic subdivision of the skeleton into its branches, permitting the distinct features of the nasal cavity topology to be individually identified.

Finding the local maxima is a sensitive process, potentially generating spurious or insignificant branches. These are mostly due to the linear connectivity of the perimeter and dense discretization of the sampling grid, and an automatic pruning method is used to eliminate spurious branches. For this, a gradient measure G=(branch length)/(skeleton radius, i.e. passage width at the base of the branch) is evaluated; branch removal occurs if G<1.1, the threshold. The threshold chosen is somewhat arbitrary and, though successful for these cases, may be different with other data. Other pruning methods are described by Weeks (1996). Finally, a small degree of Laplacian smoothing (Kobbelt et al. 1998) is applied to remove the fine ‘staircasing’ of the branch skeleton curves.

In this way, the stack of two-dimensional contours can be transformed to a stack of smooth skeleton branches with an associated thickness distribution. The process can be reversed, with closed contours again reconstructed and allowing surface mesh reconstruction via an implicit function formulation, as described in Gambaruto et al. (2008).

Figure 2b shows the right airway geometries of three subjects, models R-1, R-2 and R-3. Summary geometric characteristics are described in table 1.

The central portion of the bounding airway geometry of subject 1 has been removed, to reveal the skeleton of the airway at consecutive slices. The middle meatal skeleton branch (automatically identified as described), the upper medial passage branch and the lower meatal and medial airway branches are shown with different colouring. Separately combining the skeletonized airways and passage calibre allows an ‘averaged’ airway to be reconstructed as shown on the right. The procedure ensures separate treatment of each airway—note the averaged geometry incorporates a short superior meatus: though this is deliberately omitted from model R-1, it is present in R-2 and R-3.

The form of the airways can thus be described by the stack of the branched skeleton curves, and the techniques that characterize the skeleton curves effectively characterize the airway. For characterization, a compact representation of the skeleton branch curves is required; since they are planar, from the Frenet–Serret relations of differential geometry, they can be reconstructed given an initial tangent and curvature distribution. To illustrate, figure 2c shows the variation of curvature (κ) versus arc length (s) along the lower meatus (lm) branch of subject R-3. The mean curvature evaluated over 12 equal segments, together with the total length of the branch and the direction of the initial tangent, was found to be capable of accurately specifying any branch curve. The 14 such parameters that define the skeleton curve constitute a spatial signal, with different subject geometries being thought of as snapshots of the temporal evolution of the signal. Proper orthogonal decomposition (POD) analysis (for an introduction, see Liang et al. 2002) can then be applied to the snapshots to extract the proper orthogonal modes. Any individual branch curve can be well approximated by combining just a few (the most energetic) modes, so the coefficients of the dominant modal approximation effectively provide a characteristic trace of the geometry.

A large dataset is needed to apply the procedure usefully, but such data are not yet available given the burden of segmentation and reconstruction. As described in Gambaruto (2007), the method was instead tested using artificial data. A 243-member dataset, corresponding to a broad range of plausible lower meatus branch skeletons, was generated by combinatorial additions to a polynomial expansion of the signals defining the three subject anatomies. POD analysis was then applied and the proper orthogonal modes extracted. The coefficients for modal expansion of any curve in the dataset are found by least-squares solution; as shown on the right of figure 2c, combining just four dominant modes provides a good approximation to the lm branch of subject 3.

Complete airway characterization requires modal expansion both for all major branch skeletons and their calibre distributions; these may be treated independently or taken together to preserve feature correlation; both are straightforward in principle. Further tests and extensions to the method are under investigation, though other approaches that retain the whole airway form, based for example on Fourier descriptors, are also being evaluated. At this stage, while it can be stated that compact means to describe and characterize airway geometry exist (with skeletonization one possibility), further work is clearly needed to establish the optimal procedure.

To investigate how geometry variations affect flow, computations for the four models of figure 2b, i.e. the three replica and one averaged, were performed (Gambaruto 2007; Gambaruto et al. in press). Results indicate that both gross flow features and sensitive measures of the average geometry (respectively overall pressure loss and statistical distribution of particle residence time) lie within the range associated with the individual geometries. This cannot be taken for granted since the governing equations are nonlinear, the flow is not constrained to respond simply in proportion to changes in geometry. It is worth noting that averaging tends to regularize features—for example, the averaged airway displays a straighter and more regular septum, and so flow in the averaged airway may in fact be significantly different from that in any real subject.

Nevertheless, averaging may help establish modes of departure of the anatomical form of individuals from a population average or typical geometric modes associated with a given population (morphological variations are discussed by Churchill et al. 2004). Given a means, not only to identify such modes but also to provide a compact representation, investigating the influence of various geometric features on cavity flow can be facilitated. Furthermore, the geometry representation techniques can be applied to the structure as opposed to the cavity airspace, and provide a means to model structural alterations.

4. Internal and external anatomical form and nasal airflow

Computational modelling commonly assumes that either a blunt or fully developed velocity profile is imposed at the nares, while experimental models rely on some form of pipe inlet. In reality, the external nose influences the approaching inspired flow. Preliminary studies on the variations in inflow conditions are reported in Franke et al. (2005). In more comprehensive studies to be reported elsewhere (Taylor et al. in preparation), the geometric contraction from nasal vestibule to nasal valve is shown to reduce the dependence of the flow in the cavity to that pertaining at the nares. While the inflow condition significantly affects cavity flow, the overall flow pattern appears to be dominated by the internal nasal valve (size and orientation), together with the internal cavity geometry.

The field of view when imaging patients is generally restricted to the region of clinical interest, to limit the radiation dose. Nevertheless, to replicate the physiological situation as closely as possible and indicate the spatial origin of inspired air, the external nose and at least a local portion of the face should be modelled. While the external nose and part of the head is usually included and can be segmented from the data as discussed, it may be necessary to provide simple approximations to the portions of the face beyond the imaged region. For the case shown, these were added using CAD software (Rhino) in order to provide a free inflow to the nares. The inflow computational boundary was removed to a distance of the order of 20 cm to the front and sides of the external nose.

Simultaneous computation of the external and internal airflow is obtained by imposing different pressures at the far inflow, and at the nasopharynx, mimicking the physiological condition. Solution of the steady incompressible Navier–Stokes equations for a flow rate of 100 ml s−1 was obtained with the CFD package Fluent v. 6.3, where a third-order (MUSCL) interpolation scheme was selected to solve the momentum equation. Scaled residuals were monitored and convergence of continuity was set to 10−7, with velocity residuals typically one to two orders lower.

The paths of marker particles (neutrally buoyant and of negligible size) entering and transiting the right nose are shown in figure 3, assuming a total inhalation volume of 500 ml with an equal inspiratory rate of 100 ml s−1 for each side of the nose. The velocity contours and streamlines in a mid-sagittal plane (figure 3c) indicate rapid decay of the externally induced velocity field. Beyond the immediate vicinity of the nose, the flow resembles a simple sink, with a virtual origin located just below the naris. (The results suggest that the external boundary could be moved quite close to the naris, possibly O(5 cm) without affecting the main cavity flow. However, the computational effort to resolve the external flow is comparatively modest, as a coarse mesh suffices for the low velocities.)

Pathlines of flow, and external velocity in single inhalation. (a) Mass-less particles depicted as they transit the right nasal airway (model R-1a), at a simulated flow rate of 100 ml s−1 (through each nostril). (b) Little mixing is evident; flow entering near the tip of the nose reaches the olfactory cleft, and flow originating near the face passes quickly through the lower passages. (c) External velocity contours (m s−1) shown for a slice approximately aligned with nasal septum.

Particles arriving at the right naris are traced back to the external location from where they would originate in a single, steady inhalation. Marker particles are assigned one of three colours, depending on their position in the anterior–posterior direction at the naris (nostril entrance). The relative lack of mixing within the cavity is evident, with largely different trajectories followed by particles of different colour and hence naris origin. The different streams are not perfectly segregated, since a small recirculation in the nasal vestibule provides a limited degree of mixing. Colour coding shows that flow reaching the upper cavity predominantly passes through the anterior part of the nasal vestibule, as reported previously (Zhao et al. 2004) for a different geometry. However, figure 3 also shows that in the farfield, the upper cavity flow streams originate from directly in front of the nose, rather than from the lateral or vertical margins.

5. Mechanics of flow impact on middle turbinate

The inferior and middle turbinates (figure 1), being large and covered in erectile tissue, are key determinants of flow. The anterior portion of the inferior turbinate lies by the origin of the inspiratory flow jet that is funnelled through the nasal vestibule. It is thus involved in nasal valve function and may indeed severely limit the flow entering the cavity. The anterior head of the middle turbinate, being displaced in both posterior and superior directions, is sited to intercept the inspiratory jet. Impact of the flow on the middle turbinate generates fresh boundary layers, accompanied by high values of WSS and concomitant high levels of heat and mass exchange. This part of the nasal cavity is reportedly involved in many occupational exposure pathologies (Morgan & Montiocello 1990). The intrusion of the middle turbinate into the inspiratory jet generates complex flow patterns that have not hitherto been explored in detail. In part, this is owing to computational cost—resolving fine structures requires even finer mesh detail.

In figure 4, WSS distributions (measured in pascals) and the near-surface streamline behaviour in the region of the middle turbinate are shown for the geometries R-1 (figure 4a) and R-3 (figure 4b); mesh parameters are as outlined previously. Circumferential wall shear distributions around the perimeter of the mid-cavity slice (coloured triangle markers on surface view; cf. figure 1) are shown below the respective surface views, also with the corresponding distribution for idealized model I-1 in figure 4c. The surface contour plots are taken from a viewpoint directed upwards from the nasal valve, with the septal wall removed to reveal the middle turbinate—hence the geometry appears foreshortened. For model R-1 (figure 4a), the middle turbinate intrudes into a relatively rapid flow (max. velocity approx. 3 m s−1) due to the small nasal valve area and narrow passage calibre. A compact stagnation point is observed just below the leading edge of the turbinate, downstream of which can be seen an attachment line, as most of the impinging flow is deflected around the underside of the turbinate. High shear stress values are induced on each side of the attachment line, as fresh boundary layers develop. Bearing in mind the lower flow rate considered here, these peak shear stress values are comparable with those estimated in the computational studies of Elad et al. (2006).

Impact of inspiratory flow on middle turbinate and mid-cavity WSS. (a) (i) Lateral view of WSS distribution for model R-1. Septum removed to reveal middle turbinate; line stagnation of inspiratory jet as it impinges and flows across the turbinate visible. Fresh boundary layers induce raised shear stress—note peaks on leading edge of the intruding turbinate. (ii) WSS stress distribution along the perimeter of the mid-cavity slice (cf. figure 1); colour coding of triangles indicates location, with landmarks (a–e) also included. (iii) Comparison of scatter in predicted WSS for models R-1a,b,c; solid line represents average for all R-1 geometries. (b) (i) Model R-3, same as (a), but for subject 3. (ii) Increased patency of subject 3 geometry lowers wall shear, with peak values of approximately one-third of those for subject 1. (c) Model I-1, WSS distribution at comparable slice location for idealized model I-1 where the straighter septum has reduced the WSS levels.

Considering the circumferential wall shear variation (figure 4a(ii)), elevated values occur around the lower margins of the middle turbinate, and on the opposing and lower surfaces of the inferior turbinate and septum. Low wall shear is found in the upper reaches of each meatus.

Figure 4a(iii), representing model R-1, compares the circumferential wall shear variation for the geometries R-1a, R-1b and R-1c associated with subject 1. The slight differences in model geometry evidently result in some local scatter, as shown by point markers. However, the average distribution of wall shear (continuous line) is similar to that in figure 4a(ii). The circumferentially averaged mean WSS for R-1a,b,c at this location are, respectively, 0.0515, 0.0510 and 0.0549 Pa (overall mean 0.0524 Pa). This suggests that small differences in geometry due to subjective interpretation of image data do not affect the evaluation of a specific model geometry.

Corresponding results for a highly decongested nose are shown in figure 4b (R-3). With greater passage clearance on either side of a smaller middle turbinate and decreased jet velocity due to a larger nasal valve, wall shear levels are reduced by more than 50%. Furthermore, the surface streamlines display a more complex pattern as the impact of the inspiratory jet is more evenly distributed over the cavity walls. Regions of separation and attachment on the lateral walls are evident, and an attachment line along the lower portion of the middle turbinate is as found in R-1. The circumferential WSS displays a comparable (though reduced in magnitude) distribution to that of R-1, with low values in the olfactory cleft and in the upper parts of the meatal passages; the circumferential mean wall shear at this location is 0.0178 Pa.

Sensitivity of flow to form must be considered in attempts to reduce model complexity (e.g. Horschler et al. 2003) and perform computational modelling in a simplified nasal cavity. In this regard, computations for the idealization I-1 of model R-1 (cf. figure 2) indicated that the total nares to nasopharynx pressure drop was reduced (approx. 32%) and concomitantly the main cavity flow was significantly altered. As more than half the pressure drop occurs downstream of the nasal valve, pressure drop per se does not sensitively indicate flow conditions in the main cavity. Interestingly, for model I-1: the stagnation point for the large recirculating region of model R-1 (cf. figure 3b) was displaced anteriorly; there was a slightly lower nasal jet velocity (although nasal valve area was unaltered, the valve shape is simplified and better filled by flow); and a higher flow rate to the olfactory cleft was noted. The improved flow to the olfactory cleft is reflected in the higher wall shear there (I-1; figure 4c), compared with that found at the same location for the various approximate replica models. Elsewhere, the simplified model I-1 displays lower peak values of wall shear and lower circumferential mean (0.038 Pa).

6. Dynamics of nasal cavity flow and impact on middle turbinate

Modelling of nasal airflow generally assumes flow to be quasi-steady, though cyclic variations have been considered by Shi et al. (2006) and others. As the duration of a normal inhalation is typically an order of magnitude greater than the mean time for particles to transit the nasal cavity, quasi-steady assumptions appear reasonable. However, even where recourse is made to unsteady computations, there are a number of deficits in current modelling. At very high inspiratory rates, assuming a fixed airway boundary may not be appropriate (Brugel-Ribere et al. 2002), though detailed data are lacking. Leaving aside possible boundary deformation, current models do not attempt to resolve flow unsteadiness, potentially significant at higher inspiratory rates and in sniffing.

Detailed study of the modes of breakdown of nasal airflow has not been reported before, neither in other experimental studies nor computationally using CFD. Computations are demanding, and require careful validation; thus experimental models offer some advantages. Replica experimental models, derived from in vivo imaged data have previously been used to investigate the nasal airflow patterns and the variations in velocity (e.g. Kelly et al. (2000) and Kim & Chung (2004), among others). Here, they were applied to examine the development of flow instability at sufficiently elevated flow rates that assumptions of purely steady flow are doubtful (figure 5).

Dye visualization depicting flow irregularity and middle turbinate impact at higher inspiratory rates. Three sets of dye filament injection images depict onset of unsteadiness of inspiratory flow in three optically clear silicone phantoms. (a) (i)–(iii) Relatively stable flow in model R-3 at 150 ml s−1. The large passageway calibre reduces velocity, delaying the appearance of instability until the filaments reach the upper parts of the airway; this is shown using two sequential magnified views of the reattachment of the inspiratory jet. (b) (i)–(iii) The instability of flow in model R-1a at 170 ml s−1. Though the dye filaments appear dispersed and well mixed, flow is not truly turbulent. The magnified sequential images (ii), (iii) indicate dramatic effects on mixing and dispersion due to minor fluctuations of the inspiratory jet near the stagnation on the middle meatus. (c) (i), (ii) Filaments as they pass through the idealized model I-1 at 150 ml s−1. Shear-layer instability at the edge of the inspiratory jet is shown by successive roll-up of the dye in this region, with a detail of this region magnified in (ii). There is enhanced dispersion of dye in the anterior cavity due to the jet boundary instability and the variation in the orientation of flow approaching the middle turbinate. Note the relatively undisturbed trajectory of the major portion of flow until near the middle turbinate where filaments can be seen to oscillate as they pass through the region containing the olfactory receptors, and as they encroach on the middle turbinate.

Creation and validation of the twice life-scale replica models used are described in Doorly et al. (2008), as are procedures for dye visualization and PIV. The methodological details of most relevance are as follows. Firstly, high-resolution CT imaging showed little accumulated geometrical error accrued during fabrication, with a mean error of 0.29 mm at 1 : 1 scale. Secondly, using water as the model fluid for dye visualization, a maximum sedimentation velocity of less than 0.2 mm s−1 was obtained by adjusting the dye buoyancy, which is negligible in comparison with mean flow velocities three orders of magnitude greater. Though optical distortion occurs to varying degrees using water (ref. above), such distortions do not obscure the character of the flow revealed by dye.

As reported in Doorly et al. (2008), even at the low inspiratory flow rate of 115 ml s−1, slight spontaneous flow irregularity occurs in model R-1. Dye visualization images for a higher imposed steady inspiratory flow are compared for the two replica subject geometries R-3, R-1 and for the idealized model geometry I-1 in figure 5a–c. The images for models R-3 and I-1 correspond to a rate of 150 ml s−1; that for R-1 corresponds to 170 ml s−1. For R-1 and I-1, Reynolds number (Re) increases from approximately 900 at 100 ml s−1 to 1350 at 150 ml s−1 and to just over 1500 at 170 ml s−1; for model R-3, Re rises from 675 at 100 ml s−1 to just over 1000 at 150 ml s−1.

In all cases, flow was observed to be predominantly laminar over much, or even most of the cavity, as indicated by close scrutiny of the dye visualization images in figure 5. Considering figure 5a (model R-3), the greater nasal valve cross-sectional area lowers the inspiratory flow velocity so the flow appears to be little disturbed throughout. As evident in the pair of successive images in figure 5a(ii)(iii), the beginnings of slight irregularity are just manifest in the travelling undulation of the dye filament approaching the upper cavity. There appears to be little mixing of flow in the olfactory cleft.

For geometry R-1, a much faster jet flow is injected into the cavity at any corresponding inspiratory rate. The flow irregularity (figure 5b) is thus more noticeable, even allowing for an increased flow rate of 170 ml s−1. Patterns of dye dispersal can, however, convey misleading impressions of well-developed turbulence. Examining the flow impacting the middle turbinate (magnified details in figure 5b(ii)(iii)), it appears that rapid but small fluctuations in flow approaching the stagnation point can lead to very effective dispersal, with flow paths both diverging and bifurcating, to flow through either the middle or medial meatus. Thus although instability is evident, the form of the nasal airways appears to enhance the consequences for mixing, via augmented dispersion.

For geometry I-1, the straightened septum and the more uniform cross sections imposed during geometry simplification result in more clearly defined and visually evident instability phenomena. In figure 5c, shear-layer instability at the edge of the inspiratory jet is revealed by dye accumulation in regularly spaced patches. A magnified detail (figure 5c(ii)) indicates the dye entrained by successive rolled-up vortices in the jet edge shear layer. Also evident is the enhanced dye dispersion in the anterior cavity due to the jet boundary instability and variation in the orientation of flow approaching the middle turbinate. Note the relatively undisturbed trajectory of much of the flow until near the middle turbinate, where filaments can be seen to oscillate as they pass through the region containing the olfactory receptors, and as they encroach on the middle turbinate.

High-speed PIV measurements, and direct numerical simulation, were used to examine the flow instabilities in more detail. As described in Doorly et al. (2008), an aqueous solution of sodium benzoate was employed as working fluid instead of water for PIV, eliminating refractive index mismatch. The PIV images were acquired at a rate of 1000 fps, using a pulsed diode laser as light source. This corresponds to an equivalent rate for a 1 : 1 model in air of approximately 14 000 fps, thus capturing the transient dynamics.

Figure 6a,b is directly comparable with figure 5c, presenting both experimental (figure 6a) and computational (figure 6b) results for flow in idealized model I-1 at a nominally steady inspiratory rate of 150 ml s−1. The experimental measurements (figure 6a(i)) reveal that the time-averaged velocity in the inspiratory jet exceeds 3.5 m s−1, while weaker recirculating velocities, typically 0.5–1 m s−1, are observed in the upper anterior cavity. A pair of successive, instantaneous velocity measurements, in the region illustrated by a black rectangle, identify the unsteady nature of the flow upstream of the middle meatus. During the lapse of 0.71 ms (scaled to in vivo time) between magnified velocity images, the vortical structure at the bottom of the image has moved upwards, and regularly shed structures can be observed. The vortex spacing conforms to the patch spacing on the preceding dye visualization, and plotting the planar vorticity immediately picks out the shear-layer and the roll-up instabilities.

Steady-flow instability and sudden inspiration. (a) Experimental: (i) time-averaged PIV measurements at steady inspiratory flow rate of 150 ml s−1 through the idealized model I-1, superimposed on a particle image. The inspiratory jet exceeds speeds of 3.5 m s−1 as it is funnelled into the main cavity, with weaker recirculating velocities typically 0.5–1 m s−1 indicated by contours and direction indicators in (i). (ii) Instantaneous velocity measurements with a lapse of 0.71 ms (in vivo equivalent interval) are shown with a vorticity map derived from (i), extracted from the region illustrated with a black box in the large flowfield image. Regularly shed structures can be observed. (b) Computational: direct numerical simulation shows similar features to experiment. Instantaneous vorticity field indicates unsteadiness of flow in upper and anterior cavity, elsewhere fluctuations not evident. (c) High-speed PIV and (d) direct numerical simulation, showing the formation of the starting vortex in the idealized model as the flow is impulsed from stationary to 200 ml s−1 over a 20 ms interval. The effectiveness of mixing in the nose is greatly enhanced by transient flow dynamics, such as the sharp sniff shown.

Subsequent direct numerical simulations (lower strip) confirm the PIV and dye experimental results, where a slice through the instantaneous vorticity field again shows the unsteadiness in the jet shear layer, and the consequent irregular vorticity swept into the upper anterior part of the cavity (see again figure 5c(i)). More detailed examination of the computed velocity and vorticity fields (figure 6b(ii)) agrees closely with the PIV, bearing in mind the averaging effect associated with finite light sheet thickness is not accounted for in the computations. The jet edge instability clearly identified displays a regular frequency scaled to in vivo conditions of approximately 460 Hz. Elsewhere, the computed flow is quite steady as previously qualitatively indicated by the dye visualization. (Flow solved with Fluent v. 6.3.26, parameters: second order-in-time non-iterative advancement; second-order pressure, third-order (MUSCL) momentum approximations, respectively; hybrid mesh with four-wall layers, 0.65 million prisms, 3.6 millions tetrahedra, first layer height O(0.04 mm); amplitude 1% broadband random disturbances added to inflow; computational time step 0.01 ms.)

The dynamics of airflow and transport in sniffing is of increasing interest, particularly from a diagnostic and therapeutic viewpoint. Mild sniffing comprises short sequential bursts of increased inspiratory air flow rate (higher than 500 ml s−1), the increased pressure drop induces some degree of collapse of the external nose and probably affects the flow through the nasal valve into the main cavity. Deep sniffing produces a marked degree of alar collapse and very high transient flow rates (probably higher than 1000 ml s−1).

Recently, Zhao et al. (2004) have modelled transport and uptake during sniffing. They presented essentially quasi-steady predictions, comparing the laminar and turbulent flows, in the latter case using conventional two-equation Reynolds-averaged modelling. They found that deposition efficiency to the nasal receptor depends not only on the rate of odorant flux to the receptors, but also on its time course, with lower rates of sniff potentially increasing deposition. However, they recommended that fully unsteady calculations are attempted, particularly to resolve the Lagrangian history of inhaled species. Certainly, both the dye visualization and the PIV measurements presented here, which demonstrate irregular flow and amplification of dispersion by interaction with turbinate geometry, seriously question the validity of turbulence modelling for such flows.

What happens during rapid flow onset is of interest for olfaction and potential applications centred on inhalation. A sudden sniff, potentially augmented by some device, may induce transient flows that assist transport. The time taken to initiate significant airflow velocities in the nose can be very short; our own preliminary investigations using a hot wire probe placed below the nostril indicated flow can rise significantly well within 40 ms.

In figure 6c,d, experimental (figure 6a) and corresponding computational (figure 6b) results for the development of flow in a model nose (I-1) are shown, during a very fast (20 ms air-equivalent) ramp of inspiratory flow to a maximum of 200 ml s−1. The characteristic roll-up of the vorticity shed by the nasal valve is clearly evident. Strong transient mixing is provoked by the transit of the vortex, with the flow being more directed into the recirculation region before the vortex becomes fully established. Computations also show roll-up of the vorticity in the nasal vestibule before flow is fully established.

7. Discussion and conclusions

Procedures for modelling in vivo nasal geometry have been outlined, with attention drawn to deficits in current approaches, along with various techniques to investigate more reliably the details of form and function. Variations in model geometry due to (i) uncertainty in airway definition, (ii) idealization versus real morphology and (iii) different subject morphologies have been considered.

Uncertainty due to imperfect imaging remains an issue. Although CT technology is improving, its use is restricted and non-invasive MRI is generally less well resolved. The significance of static airway boundary uncertainty for the assessment of overall respiratory function fortunately appears to be limited from the CT-derived results presented here. However, this may not hold for more sensitive measures, e.g. transport to the olfactory airways; addressing such sensitivity is thus important. Furthermore, no effort has been made to alter images by intentionally thickening the surface to represent either the possible soft tissue displacement or the regional engorgement. Better in vivo data on dynamic airway changes are needed, but obtaining these will be challenging.

Concerning (ii) and (iii), it is evident that there are both wide inter- and intra-individual variations. Airway skeletonization was described and illustrated as one technique to decompose the geometry into essential topological connections; this enables objective reconstruction of parts of the whole architecture, and different architectures to be combined or contrasted. For brevity, other techniques based on Fourier descriptors or harmonic mapping have not been considered; the procedure used for geometric characterization is, however, secondary to its use. Future work will apply geometric deconstruction to examine modes of variation and their functional consequence.

With regard to flow, the external anatomy is often neglected. This was shown to be of importance in terms of the origin of entrained air; findings of other work were reported to indicate that consideration of the external geometry is needed to ensure modelling reflects the appropriate bias of flow to the internal nasal valve and beyond. Clearly, inflow conditions need careful consideration for bulk mixing and wall mass transfer processes, including olfaction, as opposed to mere flow resistance. Beyond this, future modelling should incorporate the interaction between the external nasal structure and the flow during deep inspiration, and during sniffing.

Complex patterns of wall shear about the middle turbinate were revealed, with very fine meshes capturing the detailed flow features around highly curved regions—particularly the front of the turbinates. WSS distributions displayed similar gross patterns in area-equivalent geometries, but comparing the straightened, simplified idealized model and the replica geometry, noticeable differences occurred in the olfactory cleft. Future work may consider mechanotransduction of stress imposed by airflow on the cavity wall.

The nature of flow unsteadiness was investigated at higher flow rates. Shear-layer instability of the jet entering the cavity was found, initiating or concomitant with temporal variation in flow pathlines near the stagnation location on the middle turbinate. The question of whether geometric simplification aids or hinders the development and amplification of flow instability is interesting and needs to be investigated. However, the effectiveness of mixing within the cavity, and regional transport, e.g. to the olfactory cleft, would appear to depend not only on particle pathline variability due to local unsteadiness, but also more subtly on the interaction between temporal and spatial variabilities in the flow. Hence, the modelling of flow and transport in the nasal cavity at flow rates where instability leads to spontaneous irregularity in pathlines using conventional pseudo-steady turbulence models is physically ill founded, and may produce spurious results.

Finally, fast inspirations or sniffing-like actions involve a clearly detectable starting vortex that originates at the expansion after the nasal valve. Its existence calls into question the reliability of assuming quasi-steady-flow conditions for fast or transient breaths and the influence of the flow on gas and particle transport processes and interactions with the walls. More work is required to assess the significance of such flow phenomena for potential applications.

Acknowledgments

We are indebted to the Biotechnology and Biological Sciences Research Council (BBSRC) for funding this research (grant nos. BBE0234441 and E18557), and to the EPSRC equipment loan pool, Oxford Lasers, the HPCx national computing resource and Fluent Ltd for additional support. We are grateful to Dr R. Almeyda for the assistance in model segmentation, ENT department, St Mary's Hospital, Paddington, and to the radiology department at that hospital for CT validation of the fabrication process used to create replica models.

Footnotes

One contribution of 11 to a Theme Issue ‘The virtual physiological human: building a framework for computational biomedicine II’.