The nasal cavity has the function of conditioning the air inhaled into the lungs by heating, humidifying, and filtering dust and virus-borne aerosols. Analyzing the flow field in the nasal cavity is vital because its function is strongly related to flow dynamics. Due to experimental limitations posed by the complex internal geometry of the nasal cavity, most previous studies have utilized Reynolds averaged Navier–Stokes based computational fluid dynamics (CFD) simulations. In this study, the flow field in a post-operative nasal cavity was evaluated using not only CFD simulations but also four-dimensional magnetic resonance velocimetry. The study was conducted under resting breathing conditions in pre- and post-operative models of a patient who received septoplasty and turbinoplasty. The experimental results confirmed balanced flow rates in the left and right nasal cavities after septoplasty and a decrease in velocity after turbinoplasty with a reduction in regions with vortices and reverse flow. Upon comparison, CFD results using the laminar, k–ω, and shear stress transport models were deemed to be consistent with the experimental results. However, there was a relatively large deviation observed with the k–ε model. Using the validated laminar CFD model, it was shown that the pressure and wall shear stress decreased after surgery.

The nasal cavity is the uppermost part of the respiratory system responsible for air conditioning (e.g., heating and humidifying) and removing fine dust or virus-contained bioaerosols. Large clusters of olfactory receptors in the nasal cavity's rear region facilitate the sense of smell. These functions are closely related to a wide range of flow phenomena such as merging and bifurcating streams, recirculating flow, and Dean-type vortices. The integrated effect of such fluid dynamic behaviors is strongly influenced by the intricacy of the nasal geometry.1 Buxton and co-workers reported that nasal ecogeographic features can vary depending on the nasal structure.2–4 In addition, Liener et al. discovered that the protrusion of the concha increases the mucosal surface area to facilitate heat and moisture exchange.5 The airflow also becomes more turbulent in the nostril as it is directed downward through a tortuous path.6 

Nasal airway obstruction interferes with smooth breathing and weakens nasal functions. In general, a deviated nasal septum (DNS) or chronic hypertrophic rhinitis (CHR) is the most common problems in rhinology. DNS is a disorder in which a curved nasal septum causes sinusitis, rhinitis, and snoring. As a typical treatment for a deviated septum, septoplasty is an operation to correct the curved area by excising the curved bone or cartilage in the nose. In CHR, the turbinates are in a chronic inflammatory state of chronic rhinitis and can be treated via turbinoplasty. After septoplasty or turbinoplasty, rhinomanometry and acoustic rhinometry are commonly employed to evaluate the therapeutic effect. Rhinomanometry measures the transnasal pressure to evaluate the nasal resistance, whereas acoustic rhinometry uses sound waves to measure the variations in the cross-sectional area. It should be noted that these two evaluation indicators are not directly correlated with nasal obstruction, and nasal function is also directly related to the flow, not the anatomical structure.7 Wiesmiller et al. investigated the effect of surgical improvement based on in vivo measurements of nasal temperature and humidity before and after septoplasty.8 However, to properly understand the mechanism of improvement in temperature and humidity after surgery, the local flow pattern must be analyzed.

Despite the complex anatomy of the nasal cavity, several scholars have attempted to experimentally measure the intranasal flow. Early in vitro studies obtained qualitative results by observing nasal air flow using smoke in air or dye in water flowing through nasal passages cast with transparent materials.9–11 Swift and Proctor measured the inspiratory flow field using a pitot tube.12 Thereafter, Girardin et al. measured the quantitative velocity field using laser doppler velocimetry.13 Schreck et al. conducted pressure measurements, hot-wire anemometry, and flow visualization using dye in water.14 Hahn et al. fabricated a right nasal model of 1:20 enlarged scale to measure the velocity profile using a hot-film anemometer probe.15 With fabrication techniques based on enhanced computerized tomography (CT) scans, detailed flow patterns have been visualized using particle image velocimetry (PIV).16–19 Most studies need to utilize refractive index matching in order for the laser sheet to pass through the complicated model. Traditional PIV obtains a two-dimensional velocity field by irradiating a laser sheet, while three-dimensional flow can be obtained via scanning stereo-PIV or tomographic PIV.20,21 However, most previous PIV research has been limited to 2D measurements, which convolutes assessment of 3D flow features such as vortices.

Alternatively, in silico computational fluid dynamics (CFD) has been actively used as an analysis tool for intranasal flow.22 Several CFD studies discussed the correlation between flow and pathology in post-operation nasal models.23–28 Ozlugedik et al.24 studied the influence of septoplasty and partial lateral turbinectomy on the flow. Although the enlarged cross-sectional area of the nasal cavity after surgery decreases the overall flow velocity due to the significantly reduced nasal resistance, the local airflow increases through the right middle meatus. Na et al.28 performed CFD analysis on an inferior turbinate surgical model to analyze the local variations in velocity, temperature, humidity, pressure gradient, and wall shear stress as well as global parameters such as pressure drop and flow rate ratio. This analysis also evaluated post-surgery flow improvement.

Although numerous computational studies have analyzed the effects of different parts of the nasal cavity and flow improvement before and after surgery, only a few studies have experimentally validated the results. Zhang and Kleinstreuer29 validated human oral airway CFD using laser Doppler velocimetry (LDV) measurement results. However, the flow geometry was assumed to be a simple axisymmetric constricting tube, and LDV also has a limitation in that it is a point measurement. Bass and Lingest30 investigated the deposition of micro-particles in the upper airways and experimentally validated a two-equation turbulence model. However, they also assumed the flow geometry to be a simple 90° circular cross section with bends, and the experimental data used for CFD validation were aerosol deposition data, not flow field data. Chen et al.31,32 calculated the deposition, evaporation, and growth of aerosols according to flow rate, temperature, relative humidity, and wet-dry airway wall conditions in the human airways and experimentally validated CFD. However, they also assumed the flow geometry as a simple constricting bend, and the data used for CFD validation were aerosol deposition and droplet growth. They were not able to directly compare the flow field.

Magnetic resonance velocimetry (MRV) utilizes phase-contrast magnetic resonance imaging (MRI) and can measure the flow field of three-directional velocity components in three dimensions. MRV is useful for studying various vascular flow and engineering applications, because it can measure the flow velocity without optical access.33,34 Four-dimensional (4D) MRV enables the evaluation of the temporal variations of a 3D flow field in a given period.35 Therefore, this technique has been applied to biofluid measurements within vascular and respiratory systems.36–38 Banko et al.37 examined the flow through the larynx, trachea, and bronchus under oral breathing conditions. Jalal et al.39 used 4D MRV to investigate the velocity field from the middle of the trachea to the bronchus. Collier et al.40 utilized MRV as a validation tool for CFD simulations of the 3D mean flow in the central airway tree. Borup et al.41 used MRV to measure the mean flow field within the nasal cavity. Nonetheless, the dearth of existing research using 4D MRV for nasal measurements renders this topic worthwhile to study.

The primary objective of this study is to examine the feasibility of 4D MRV for nasal flow research by investigating the complex 3D flow structure to quantitatively assess the results of nasal surgeries. The cyclic nasal flow structure before and after combined septoplasty and turbinoplasty is elucidated. The local 3D flow patterns at all temporal phases of the cycle are obtained through 4D MRV experiments and are used to validate CFD simulations. The wall shear stress and pressure, however, which cannot be obtained through MRV, are calculated from the validated CFD results. These metrics allow assessment of the improvement in flow and pressure drop after the operation. The manuscript is organized as follows: Sec. II provides an overview of the experimental methods detailing the apparatus used, and the numerical analysis undertaken. Section III presents the results, sequentially ordered based on the experiments and their corresponding numerical analysis. Finally, Sec. IV offers a discussion on the relevance of the findings and concludes with a comprehensive summary.

The patient was a 52-year-old Asian male of 180 cm height and 88 kg weight, diagnosed with a deviated septum and chronic rhinitis. The patient underwent turbinoplasty, which involved trimming the bilateral inferior and middle turbinate, and septoplasty to correct a deviated nasal septum that leaned to the left. Both operations were performed simultaneously. The septum had a spur (i.e., a bony or cartilaginous projection) that was removed. Following these procedures, the septum was straightened, and the space between the turbinates expanded, improving the airway passage. Patient-specific geometric information was obtained through CT scans for the pre- and post-operative cases. Before the operation, a third-generation multislice helical CT scanner (Lightspeed Ultra, GE Medical Systems) was used to perform noncontrast scanning at 120 kVp with a slice thickness of 0.25 mm. After the operation, a third-generation dual source CT scanner (SOMATOM Force, Siemens Healthineers) was used to perform noncontrast scanning at 100 kVp with a slice thickness of 1.0 mm. The CT images were saved in Digital Imaging and Communication in Medicine (DICOM) format. The DICOM file was segmented using Mimics Innovation Suite software (Materialise NV, Belgium) and converted to stereolithography (STL) format by 3-matic software (Materialise NV, Belgium). The study was approved by the Institutional Review Board (IRB) of Seoul National University Hospital (IRB No.2005–092-1123). Figure 1 shows the full 3 D structure used to construct the flow path. In the figure, the top and bottom images correspond to before and after surgery, respectively, and the four consecutive images on the right show the cross sections at 20 mm intervals. After surgery, it can be qualitatively observed that the area is widened and the channel is straightened.

FIG. 1.

CT scan results for (a) pre-operative and (b) post-operative model. Cross sections in the coronal direction at 20 mm intervals are shown in the four consecutive images on the right.

FIG. 1.

CT scan results for (a) pre-operative and (b) post-operative model. Cross sections in the coronal direction at 20 mm intervals are shown in the four consecutive images on the right.

Close modal

The reconstructed full-scale test section is illustrated in Fig. 2. The nasal cavity model was constructed using 3D Builder.42 In MRV measurements, eliminating small air bubbles in the test section is important for preventing the detection of undesired signals and to obtain the proper flow structure. Thus, the nasal cavity model was 3D printed using stereolithography (SLA) with a transparent material (Accura ClearVue, 3D SYSTEMS). The model is segmented into three sections (i.e., front nasal cavity, rear nasal cavity, and trachea) fastened with flanges to avoid the use of inner supports when printing overhanging structures. To ensure accurate alignment, dowel pins were installed at flange interfaces. Due to limitations in the CT scan data, facile features such as the external nose shape, which can potentially affect the internal flow, were not considered in this study.43 An ideal uniform flow was utilized for the inlet boundary condition of the nostril, for data reproducibility and comparability. Flow-conditioning sections were placed both in front of the nostril and downstream of the trachea to create uniform flow, attenuating the secondary flow caused by the curved inlet tubing and connecting components. The upstream flow-conditioning section was connected to a reservoir, and the downstream section was connected to a custom pulsatile flow pump, which was fabricated using a syringe and translation stage controlled via Arduino. The test section with connected acrylic pipes and reservoir was placed in the scan room, and the electric devices including the pulsatile flow pump and flowmeter were placed in the control room [Fig. 2(b)]. In particular, a dilute aqueous solution of copper sulfate at concentration of 0.06 mol/l was used as the working fluid. Copper sulfate increases the signal-to-noise ratio (SNR) without affecting the physical properties of water.

FIG. 2.

Experimental setup. (a) Picture of test section with flow-conditioning pipes. (b) Schematic of test section and reservoir placed in the MRI scan room, and the remaining flow-driving system including flowmeter placed in the control room.

FIG. 2.

Experimental setup. (a) Picture of test section with flow-conditioning pipes. (b) Schematic of test section and reservoir placed in the MRI scan room, and the remaining flow-driving system including flowmeter placed in the control room.

Close modal
The flow was generated to resemble normal breathing conditions at rest in air. The corresponding respiratory minute volume was calculated to be 6 l/min at 0.2 Hz respiratory rate (RR) and 0.5 l of tidal volume (TV), which resulted in a mean streamwise velocity of 0.32 m/s in the trachea at peak respiration.17,25,44–46 Dynamic similarity between Reynolds and Womersley numbers was achieved using the copper sulfate solution as a surrogate fluid replicating air. Reynolds and Womersley numbers are defined as follows:
R e = Q D H ν A ,
(1)
W o = D H 2 ω ν ,
(2)
where Q denotes the volumetric flow rate, DH indicates the hydraulic diameter, A represents the cross-sectional area of the tracheal airway, ν denotes the kinematic viscosity of the fluid, and ω indicates the frequency of oscillation. The flow conditions in water and air are summarized in Table I. A pulsatile pump was used to generate a sinusoidal waveform for the volumetric flow rate, which resembles a human respiratory cycle. The Reynolds number at the peak respiration point was 591, and the Womersley number of the cycle was 1.1 based on the mean tracheal hydraulic diameter. In this regime, quasistatic convective flow prevailed with negligible effect of oscillations.47 The inlet flow rate was confirmed using a turbine flowmeter (Tecflow, IR-Opflow) (Fig. 3). Near phase 0° and 200°, the flow rate does not change smoothly due to limitations in the piston pump hardware.
TABLE I.

Flow characteristics.

Reynolds number 591
Womersley number 1.1
Air  Respiratory rate (RR)  0.2 (Hz) 
Tidal volume (TV)  0.5 (l) 
Respiratory minute volume (MV = TV × RR)  6 (l/min) 
Mean streamwise velocity in trachea at peak respiration  0.32 (m/s) 
Hydraulic diameter in trachea (DH 14.3 (mm) 
Kinematic viscosity  1.48 × 10–5 (m2/s) 
Water  Mean streamwise velocity in trachea at peak respiration  0.02 (m/s) 
Respiratory rate  0.025 (Hz) 
Kinematic viscosity  0.96 × 10–6 (m2/s) 
Peak flow rate  0.4 (l/min) 
Reynolds number 591
Womersley number 1.1
Air  Respiratory rate (RR)  0.2 (Hz) 
Tidal volume (TV)  0.5 (l) 
Respiratory minute volume (MV = TV × RR)  6 (l/min) 
Mean streamwise velocity in trachea at peak respiration  0.32 (m/s) 
Hydraulic diameter in trachea (DH 14.3 (mm) 
Kinematic viscosity  1.48 × 10–5 (m2/s) 
Water  Mean streamwise velocity in trachea at peak respiration  0.02 (m/s) 
Respiratory rate  0.025 (Hz) 
Kinematic viscosity  0.96 × 10–6 (m2/s) 
Peak flow rate  0.4 (l/min) 
FIG. 3.

Flow rate measurement from flowmeter.

FIG. 3.

Flow rate measurement from flowmeter.

Close modal

In this study, the MRV experiment was performed using the Siemens 3T MRI system (MAGNETOM Tim Trio, Siemens Healthcare, Erlangen, Germany) installed at Seoul National University Hospital, Korea. The scan information is detailed in Table II. Each scan measured a single field-of-view (FOV) of 80 × 60 × 104 mm3 in the x-, y-, and z-directions, using a voxel size of 0.5 × 0.5 × 1.0 mm3. Slices were normal to the z-direction, and the 1.0 mm thickness was interpolated down to 0.5 mm. The final processed data, thus, had an isotropic voxel size of 0.5 mm.

TABLE II.

Scanner details.

Scanner type  3.0T Siemens Magnetom Tim Trio 
4D flow sequence  Siemens WIP sequence 
Scanner software version  VB17A SW version 
Gradients  45 (mT/m) 
Max slew rate  200 [mT/(m ms)] 
RF coil  6-channel body
8-Channel Spine 
Orientation  Coronal 
Velocity encoding (VENC) – x, y, z  [Pre-operative] 10, 20, 20 (cm/s)
[Post-operative] 8, 13, 13 (cm/s) 
Echo time, repetition time (TE, TR)  [Pre-operative] 5.88, 11.3 (ms)
[Post-operative] 6.69, 12.1 (ms) 
K-space segments  [Pre-operative] 43
[Post-operative] 40 
Flip angle  15° 
Acquisition window  39 500 (ms) 
Number of phases  20 
Single acquisition time  193 (min) 
Number of scans  [Flow on] 2
[flow off] 1 
Scanner type  3.0T Siemens Magnetom Tim Trio 
4D flow sequence  Siemens WIP sequence 
Scanner software version  VB17A SW version 
Gradients  45 (mT/m) 
Max slew rate  200 [mT/(m ms)] 
RF coil  6-channel body
8-Channel Spine 
Orientation  Coronal 
Velocity encoding (VENC) – x, y, z  [Pre-operative] 10, 20, 20 (cm/s)
[Post-operative] 8, 13, 13 (cm/s) 
Echo time, repetition time (TE, TR)  [Pre-operative] 5.88, 11.3 (ms)
[Post-operative] 6.69, 12.1 (ms) 
K-space segments  [Pre-operative] 43
[Post-operative] 40 
Flip angle  15° 
Acquisition window  39 500 (ms) 
Number of phases  20 
Single acquisition time  193 (min) 
Number of scans  [Flow on] 2
[flow off] 1 

We used a 4D flow sequence comprising 3D Cartesian spatial encoding and acquisition along with a four-point velocity encoding (VENC) for three separate VENC values in each direction. The x-, y-, and z-direction VENC values (which can be considered the dynamic range of velocity) were set as 20, 10, and 20 cm/s in the pre-operative model and 13, 8, and 13 cm/s in the post-operative model, respectively. Prospective gating was applied for the 4 D MRV, and a TTL gating signal was transmitted by the Arduino at intervals of 40 s, which also controlled the pulsatile flow pump. The acquisition window was set to 39.5 s and equally segmented into 20 temporal phases for a temporal resolution of 2.0 s. In this sequence, the repetition time (TR) was set to fit the relation 4 × TR × (number of segments) into the temporal resolution. For instance, 4 × 11.3 ms × 43 yields 1.94 s in the case of the pre-operative model. The total acquisition time was calculated as (period of one cycle) × (number of k-space lines in phase encoding direction) × (number of slices) / (number of segments), which resulted in 193 min in the case of the pre-operative model. Acceleration methods such as SENSE or GRAPPA were not applied in this study.

The open-source MATLAB code “mapVBVD” (developed by Philipp Ehses) was used to reconstruct the k-space data, which encodes the Siemens raw data for further processing in MATLAB. The data underwent a coil combination procedure in MATLAB. During this process, signals from two arrays in both the torso and spine coils were integrated using the weighted sum method for each time frame. This approach adjusts the signal values based on the coil sensitivity at each voxel. In the resulting combined data, the phase differences for each voxel are determined, and VENC values are multiplied for each direction.

Two flow-on scans were averaged, and one intermediate flow-off scan was subtracted to obtain the final velocity data. All three scans were acquired at identical VENC and protocol settings. The flow region was extracted using a volumetric mask generated from the signal intensity data in the flow-off scan, forcing the data in the outer regions as well as the wall to be zero. For each temporal phase and velocity component, the uncertainty involved in the velocity measurement was calculated from the raw data using the equation proposed by Bruschewski et al.,
σ u = 1 / 2 NSA ( Var { u 1 ( ROI ) u 2 ( ROI ) } ) ,
(3)
where NSA denotes the number of signal averages, and the factor of 2 corresponds to two datasets.48 The variance of velocity was obtained from the difference between two independent datasets containing voxels within the region of interest (ROI). The uncertainty for the different temporal phases varied as 0.12–0.14, 0.19–0.25, and 0.20–0.33 cm/s in the x-, y-, and z-directions, respectively. Overall, the evaluated uncertainty was higher in the direction with higher velocity values because of the large phase differences of the hydrogen proton spins within the voxel. Nevertheless, these values were within 2.5% of the VENC for all directions and phases.

The numerical analysis conducted in this study involves two steps. First, the CFD model is validated with the MRV experimental results by comparing the results of the incompressible steady flow simulation using water as the working fluid at a flow rate of 0.4 l/min. The velocity field, root mean square (RMS) value of the velocity components, and integral parameters were analyzed for comparison. Second, the CFD simulations are performed with air as the working fluid at a flow rate of 6 l/min using the numerical model and grid validated in the first step. Factors, such as variations in the pressure field, wall shear stress, and turbulence intensity, which could not be evaluated through MRV, were numerically analyzed. In both cases, the Reynolds number in the trachea was 591, and a no-slip boundary condition was applied to the wall.

Raw surface data from a CT scan (.stl format) was converted into a solid model using Ansys SpaceClaim. The volumetric error between the solid model and raw data was less than 0.8%. A hybrid mesh comprising a tetrahedron cell mesh and prism layer mesh near the wall was employed. To obtain a grid-independent solution, 12.8 × 106 meshes in the pre-operative model and 11.5 × 106 meshes in the post-operative model were used, as illustrated in Fig. 4. The solid model for the CFD has exactly the same geometry as that of the experiment, and the inlet and outlet are located where uniform flow is formed after flow-conditioning sections in the experimental model [Fig. 4(a)]. Steady CFD simulations were performed for flow conditions corresponding to the peaks of inhalation and exhalation in the experiment, since the oscillatory flow can be assumed quasistatic (as described in Sec. II B). Ansys Fluent was utilized on an Intel® Xeon® Gold 6348 Processor (96 cores) for parallel computation.

FIG. 4.

(a) Flow field model. (b) Flow field grid of dashed box in the left figure. (c) Coronal cross section of pre-operative flow field grid for A–A′ plane. (d) Coronal cross section of post-operative flow field grid for A–A′ plane.

FIG. 4.

(a) Flow field model. (b) Flow field grid of dashed box in the left figure. (c) Coronal cross section of pre-operative flow field grid for A–A′ plane. (d) Coronal cross section of post-operative flow field grid for A–A′ plane.

Close modal

In general, the breathing condition at rest follows laminar flow. However, as the patient was diagnosed with DNS and CHR, the cross-sectional flow area within the nasal cavity changed significantly and also asymmetrically. Despite the low Reynolds number condition, local occurrence of partial turbulence is possible.44,49 Mylavarapu et al. proposed that the k–ω turbulence model is the most suitable when comparing the simulation and pressure measurement results at relatively high flow rates.50 However, experimental validation using flow data are insufficient and the criteria when to choose certain models are unclear. Thus, this study investigated various turbulence models such as k–ε, k–ω, and SST, and also examined the laminar model as well.

The change in the cross-sectional area of the nasal cavity is an important evaluation factor for turbinoplasty. The 3D structure of the nasal cavity can be visualized based on the magnitude of the MRI signal. The variation in the cross-sectional area from the nostril to the pharynx is presented in Fig. 5. Setting the cross section perpendicular to the dominant local flow can provide a clear representation of the flow characteristics within each plane. However, this approach may lead to overlapping cross sections. As a result, we have also considered physiological features in the selection of cross sections. The cross sections were rotated from 0° (horizontal plane) to 90° (vertical plane) in the nasal vestibule region, and then, multiple 90° planes were evaluated in the turbinate region. Thereafter, the cross sections were rotated from 90° to 180° (horizontal plane) in the nasopharynx region. The pre-operation and post-operation 3D structures of the nasal cavity are depicted in Figs. 5(a) and 5(b), respectively. The cylinder before the nostril is a part of the inlet channel. It can be observed that the flow channel has been expanded after operation of the nasal cavity. The quantitative variations in the cross-sectional area are detailed in Fig. 5(c). The area in the turbinate region increased remarkably after turbinoplasty, but those of the nasal vestibule and nasopharynx region did not display any notable changes.

FIG. 5.

Evaluation of cross-sectional area using MRV data. Three-dimensional structure of nasal cavity in (a) pre-operation scenario and (b) post-operation scenario. (c) Variation of cross-sectional area from nostril to pharynx region.

FIG. 5.

Evaluation of cross-sectional area using MRV data. Three-dimensional structure of nasal cavity in (a) pre-operation scenario and (b) post-operation scenario. (c) Variation of cross-sectional area from nostril to pharynx region.

Close modal

The bold solid lines in Fig. 6 illustrate the variations in the cross-sectional area ratio between the left and right cavities owing to septoplasty, wherein the slices were sequenced according to Fig. 5. The area ratio was evaluated only in the vestibule and turbinate region. The thin dashed lines are cross-sectional areas of the left and right airways. The green lines are based on the pre-operative data, whereas the blue lines indicate the post-operation results. This graph confirms that there is individual variation in the nasal flow path. As can be seen from the thick lines, the overall symmetry of the left and right nasal cavities has been improved. Although the cross sections are not directly indicative of air flow dynamics, it provides a basis for the flow fields described later.

FIG. 6.

Evaluation of left and right cross-sectional area using MRV data.

FIG. 6.

Evaluation of left and right cross-sectional area using MRV data.

Close modal

A significant difference can be observed in the rear region of the nasal vestibule and the front region of the turbinate, suggesting that there has been some appreciable improvement in distributing inhaled air evenly along the two passages as it enters the turbinate. Moving further into the turbinate, such trend is no longer discernible, demonstrating that septoplasty has been performed in the rear region of the nasal vestibule and the front region of the turbinate. Turbinectomy increases the cross-sectional area of the turbinate, which results in reduced nasal resistance and smoother breathing.

The transient flow rate calculated from the MRV data is plotted in Fig. 7. For flow calculation, the velocity and area were obtained from the cylindrical flow channel located before the nostril. The waveform of the flow rate evaluated from the MRV data was very similar to that measured by the flowmeter (Fig. 3). The cycle is divided into 20 temporal phases, wherein phases 1–10 and 11–20 correspond to inhalation and exhalation, respectively. The flow is further delineated into accelerating, peak, and decelerating phases of inhalation and exhalation.

FIG. 7.

Flow rate calculated from MRV data. I and E denote inhalation and exhalation, and A, P, and D indicate acceleration, peak, and deceleration, respectively.

FIG. 7.

Flow rate calculated from MRV data. I and E denote inhalation and exhalation, and A, P, and D indicate acceleration, peak, and deceleration, respectively.

Close modal

The velocity fields in the 3D nasal cavity phantom are displayed in Fig. 8, with the pre-operation and post-operation models presented in Figs. 8(a) and 8(b), respectively. From top–left to bottom–right, each figure corresponds to the plane-normal velocity field of IA, IP, ID, EA, EP, and ED. The black arrow with the small square shows the positive normal orientation of the cross section. In the case of the pre-operation model, high velocity magnitude was observed in the right nasal cavity for all temporal phases. Notably, the high velocity on the right was maintained up to the nasopharynx region at the inhalation peak. In addition, the velocity in the middle and superior meatus was higher than that in the inferior meatus, contrary to a previous study which reported that the flow is predominantly distributed in the inferior and middle meatus.51 In the lower portion of the vestibule region, recirculation flow was observed during inhalation as flow separation occurred along the 90° curved passage and nostril sill. Moreover, owing to the abrupt area expansion, reverse flow existed in the transition region from the turbinate to nasopharynx (as illustrated by the change in color). During exhalation, high velocity magnitude was again observed in the middle and superior meatus. Unlike inhalation, slight flow separation occurred from the nostril in the outward direction (marked with a red dotted line). In the nasopharynx region, reverse flow is caused by the shape of the curved pipe and area expansion.

FIG. 8.

MRV plane-normal velocity fields in the nasal cavity phantom. (a) Pre-operation model. (b) Post-operation model.

FIG. 8.

MRV plane-normal velocity fields in the nasal cavity phantom. (a) Pre-operation model. (b) Post-operation model.

Close modal

In the case of the post-operation model, the velocity was more uniformly distributed on each side due to septoplasty. Compared to the pre-operation model, the maximum velocity magnitude decreased in all temporal phases because of the enlarged nasal area after turbinoplasty. Notably, the velocity distribution in the vestibule region became more uniform and similar in both the left and right nasal cavity. At the inhalation and exhalation peaks in the post-operation model, the high velocity magnitude region was displaced leftward for the middle and inferior meatus. The reverse flow in the vestibule and nasopharynx regions also decreased due to the area expansion after turbinoplasty.

The out-of-plane velocity contours with in-plane vectors corresponding to the slices illustrated in Fig. 8 (from the nostril, past the nasal vestibule region, and to the nasopharynx region) are presented in Fig. 9. The two left slices and two right slices of each sub figure correspond to the IP and EP, respectively. As observed in the nostril slice of Fig. 9(a), the velocity in the right side was higher than that in the left and nonuniformly distributed in the pre-operation case in both IP and EP. In EP, the asymmetry is quite prominent, along with reverse flow in the left side. Also, the red horizontal bar-shaped flow observed in the left nostril is a jet-like flow due to the downward flow from the upper meatus and the upward flow from the lower meatus. This possibly causes large pressure drop, which would contribute to discomfort in breathing. A vortex also occurs in the anterior region (at the top of the image) that exists even after surgery. Although the flow becomes relatively uniform after surgery, higher out-of-plane velocity is still observed in the right side.

FIG. 9.

Velocity contour with vectors. (a) Sagittal plane in nostril. (b) 90th slice, between nasal vestibule and turbinate region. (c) 120th slice, anterior portion of turbinate region. (d) 150th slice, middle portion of turbinate region. (e) 210th slice, nasopharynx region. Each slice is situated at a distance of 1.5 cm from 90th to 210th slice.

FIG. 9.

Velocity contour with vectors. (a) Sagittal plane in nostril. (b) 90th slice, between nasal vestibule and turbinate region. (c) 120th slice, anterior portion of turbinate region. (d) 150th slice, middle portion of turbinate region. (e) 210th slice, nasopharynx region. Each slice is situated at a distance of 1.5 cm from 90th to 210th slice.

Close modal

In Fig. 9(b), slight reverse flow occurred in the top-right region in IP after surgery (marked with a red dashed line), which adversely affects inhalation in the downstream anterior turbinate region. After surgical intervention, the flow field in EP displayed greater improvement than that in the IP. Notably, the reverse flow disappeared, and the left and right velocities exhibited relatively higher uniformity.

The most anterior slice in the turbinate region is presented in Fig. 9(c). The flow primarily passes through the entrance of the middle meatus and then into the superior meatus, particularly in the right cavity. Slight reverse flow can be observed in the superior meatus of the post-operation model during IP, extending from the nasal vestibule region. Notably, a vortex occurred at the entrance of the middle and superior meatus (marked with a red dashed line). It has been shown that the local shear stress increases near a vortex boundary, which causes discomfort and inflammation.52 After surgery, the upward velocity increased in IP, indicating a potential enhancement in the patient's sense of smell owing to the increased flow in the olfactory region. Although the specific location of the olfactory sensing organ varies among individuals, it is generally located at the center of the turbinate region.28 The in-plane secondary velocity in the middle and superior meatus region tends to increase during IP and decrease during EP. The increased in-plane velocity likely increases the flow rate into the maxillary sinus connected to the upper portion of the middle meatus during inhalation. Similar to the case presented in Fig. 9(b), the flow field in EP improved considerably with surgery.

As observed from Fig. 9(d), the out-of-plane velocity at the entrance of the middle meatus decreased after surgery during both inhalation and exhalation. The velocity magnitude became uniform in the left and right cavities because of septoplasty, and the maximum velocity tended to decrease after turbinoplasty. However, the in-plane velocity did not significantly vary between the pre- and post-operation cases.

The most anterior slice (#210) in the nasopharynx region is presented in Fig. 9(e). The flow became more uniform with surgery in this region even though septoplasty only altered the cross-sectional area up to about the 120th slice (i.e., anterior portion of the turbinate region) as shown in Fig. 5. After operation, reverse flow regions diminished due to area expansion, which further reduced both the out-of-plane and in-plane velocities in the nasopharynx region. In particular, the vortex in EP decreased significantly.

Several previous studies have attempted at experimentally validating CFD data of nasal cavity flow, but they either compared only pressure drop without velocity field information,50,53 or compared the flow field obtained using an invasive method such as a hot-wire anemometer probe.54 Although Doorly et al.55,56 used noninvasive PIV, it was limited to 2D measurements. This study utilized noninvasive MRV to measure the full three-component velocity field within the entire 3D nasal cavity domain. Therefore, more thorough validation could be performed for the CFD data.

The CFD results were both qualitatively and quantitatively compared with the MRV observations of the flow field. Qualitative comparison involved out-of-plane velocity contours and in-plane vectors. The flow fields for the IP and EP phases were compared in four representative regions: nostril, nasal valve, turbinate, and choana (Figs. 10–13, respectively). The experiment and CFD are in good agreement with the complex flow phenomena including non-uniform flow velocity distribution, reverse flow, secondary flow, and vortices in the pre-operative nostril (Fig. 10) and choana (Fig. 13). In Fig. 11, it can be confirmed that reverse flow measured in the lower part during inspiration and expiration before operation and in the right nasal cavity during inspiration after operation are captured in all the CFD models. In addition, for the turbinate area of Fig. 12, the velocity magnitude is higher in the middle and superior meatus than in the inferior meatus on the right side of the pre-operative model for the CFD as well as the experiment, and the high velocity regions move to the middle and inferior meatus on the left side in the post-operative model. The pre-operative and post-operative results for different RANS models are also depicted. The results using the different models exhibited highly similar trends in the normal velocity distribution, in-plane flow direction, vortex position, and reverse flow position. This is because the Reynolds stress term has little effect due to the low Reynolds number. Only the k–ε model had some slight discrepancies. In this study, the k–ε model, which has weaknesses in the boundary layer and separation flow, is not appropriate for nasal cavity analysis because Re # is small and the wall has a large influence on the flow field due to the complex shape of the nasal cavity. A specific region showing the difference due to the k–ε model will be described later (Fig. 15).

FIG. 10.

Contour for out-of-plane velocity with vectors for in-plane velocity in the nostril, which corresponds to Fig. 9(a). Size and interval of the vectors are adjusted for comparison.

FIG. 10.

Contour for out-of-plane velocity with vectors for in-plane velocity in the nostril, which corresponds to Fig. 9(a). Size and interval of the vectors are adjusted for comparison.

Close modal
FIG. 11.

Contour for out-of-plane velocity with vectors for in-plane velocity in the nasal valve, which corresponds to Fig. 9(b). Size and interval of the vectors are adjusted for comparison.

FIG. 11.

Contour for out-of-plane velocity with vectors for in-plane velocity in the nasal valve, which corresponds to Fig. 9(b). Size and interval of the vectors are adjusted for comparison.

Close modal
FIG. 12.

Contour for out-of-plane velocity with vectors for in-plane velocity in the turbinate, which corresponds to Fig. 9(d). Size and interval of the vectors are adjusted for comparison.

FIG. 12.

Contour for out-of-plane velocity with vectors for in-plane velocity in the turbinate, which corresponds to Fig. 9(d). Size and interval of the vectors are adjusted for comparison.

Close modal
FIG. 13.

Contour for out-of-plane velocity with vectors for in-plane velocity in the choana, which corresponds to Fig. 9(e). Size and interval of the vectors are adjusted for comparison.

FIG. 13.

Contour for out-of-plane velocity with vectors for in-plane velocity in the choana, which corresponds to Fig. 9(e). Size and interval of the vectors are adjusted for comparison.

Close modal
An integral parameter was introduced for more quantitative comparison. The characteristics of the streamwise and transverse velocities of the flow field in a given cross section can be described by the parameter E, defined as follows:
E = A u u · n ̂ n ̂ 2 d A A u · n ̂ 2 d A ,
(4)
where u denotes the velocity vector, n ̂ represents a unit vector normal to the cross section, and A indicates the cross-sectional area.37,57 The numerator and denominator represent the strength of the in-plane and out-of-plane velocity, respectively. Therefore, the parameter E indicates the strength of the in-plane velocity compared to the out-of-plane velocity in a given cross section. A high E indicates that the flow elements were more dispersed in the transverse direction than in the streamwise direction. The results of E parameters evaluated for the experiments and CFD simulations before and after surgery are presented in Fig. 14, with the colored regions representing experimental uncertainty. The numbers of the cross section slices on the x-axis follow the nasal vestibule region (slice #1–#90), turbinate region (slices #91–#211), and nasopharynx region (slice #212–#301) in Sec. III A. The E parameters for pre-operation IP, post-operation IP, pre-operation EP, and post-operation EP are presented in Figs. 14(a)–14(d), respectively.
FIG. 14.

E parameter from nostril to nasopharynx, with colored regions representing experimental uncertainty. (a) Pre-operation IP, (b) post-operation IP, (c) Pre-operation EP, and (d) post-operation EP.

FIG. 14.

E parameter from nostril to nasopharynx, with colored regions representing experimental uncertainty. (a) Pre-operation IP, (b) post-operation IP, (c) Pre-operation EP, and (d) post-operation EP.

Close modal

As observed from Fig. 14, the decrease in E parameters for pre- and post-operation is noticeable in the nasal valve (#80–100) and choana (#212–240). Because of the high flow resistance in the turbinate region during pre-operative inhalation, the airflow does not easily enter the turbinate region from the nasal valve. As a result, the secondary flow intensifies, resulting in a large E parameter. In the choana, as the air passes through the turbinate region with a narrow cross section, relatively higher velocity creates a stronger shear layer resulting in substantial mixing. This also results in intensified secondary flow and a large E parameter. However, since the flow resistance of the turbinate region is reduced in post-operative inhalation, air flow from the nasal valve to the turbinate region is relatively easier compared to pre-operative inhalation. Thus, the secondary flow is relatively small, resulting in a decrease in the E parameter. Assuming that the respiratory flow rate before and after surgery is the same when the air flows into the choana, as it passes through the turbinate region with a wider cross-sectional area after surgery, the E parameter is reduced as slower flow results in less mixing compared to pre-operative inhalation. For choana inhalation depicted in Fig. 13, the blue colored out-of-plane velocity contour shows that surgery reduced the area of reverse flow (depicted by yellow color), while the velocity distribution became more uniform. The in-plane vectors illustrate that the right vortex has been almost eliminated.

Similarly, in pre-operative exhalation, when the airflow enters the turbinate region from the choana, the flow resistance is large, resulting in intensified secondary flow and a large E parameter. Since the airflow entering the nasal valve through the narrow turbinate region undergoes a rapid cross-sectional area change, the E parameter increases again as the secondary flow increases. On the other hand, the left and right flow distribution becomes uniform and the cross-sectional area becomes wider; thus, the cross-sectional area change is relatively small compared to pre-operative exhalation in the choana and nasal valves. This reduces the occurrence of secondary flow and the value of the E parameter. In the red colored contours of nasal valve exhalation in Fig. 11, the out-of-plane velocity is more uniform and the in-plane vortex is reduced after surgery. In other words, assuming that the same amount of flow is used for respiration before and after surgery, the ratio of in-plane velocity to streamwise velocity decreased because the airflow traveled through a wider passage with lower resistance after surgery. This is in alignment with the more uniform post-operation flow in Figs. 10–13.

The flow in the trachea is laminar with a low Reynolds number of roughly 600, and thus, one would expect the laminar CFD model to work reasonably well in this region. However, there are various regions in the nasal cavity with high local Reynolds number arising from complex flow obstacles and rapid variations in the cross-sectional area. Thus, this study employed different turbulence models and showed that the E parameters of all four CFD models were reasonably consistent with the experimental results within experimental uncertainty. Across all measurement regions, the kω and SST models delivered essentially the same results as those of the laminar model, demonstrating that the turbulent Reynolds stress term has a minor effect on the calculations. The discrepancy at small and large slice numbers is likely due to the relatively low signal to noise ratio (SNR) at the edge of the MRI imaging domain.

In the choana [Fig. 14(a)], the E values of the kε model were smaller than those observed in the experiment and other CFD models. The reason for the discrepancy of this pre-operative inhalation case can be explained with Fig. 15, which shows a region at the end of the turbinate [i.e., choana region in Fig. 15(a)]. The z-velocity component contour and in-plane vectors of the experimental data in Fig. 15(b) display a large separation zone due to flow exiting between the inferior and middle turbinate. The laminar CFD model in Fig. 15(c) shows a similar flow structure as that of the MRV data. However, the flow structure for the k–ε model in Fig. 15(d) displayed a different trend. The flow follows the curvature of the surface more closely compared to the experiment and other CFD models. This is due to delayed flow separation. The k–ε model used the standard wall function, whereas the remaining CFD models used the scalable wall function. In general, the k–ε model cannot predict the flow near the wall as accurately as the other models. In-plane recirculation flow and vortices (in the red dashed circle) are also not fully captured. Therefore, the E parameter of the kε model in the choana [Fig. 14(a)] is smaller than that of the experiment and other CFD models.

FIG. 15.

z-direction velocity component contour and in-plane vectors at posterior portion of turbinate region for the pre-operative IP condition. (a) Observation plane from MRV magnitude data. (b) MRV experiment. (c) Laminar CFD model. (d) k–ε model.

FIG. 15.

z-direction velocity component contour and in-plane vectors at posterior portion of turbinate region for the pre-operative IP condition. (a) Observation plane from MRV magnitude data. (b) MRV experiment. (c) Laminar CFD model. (d) k–ε model.

Close modal

Assuming that the respiratory flow rate is the same in the resting state before and after surgery, the variations in the respiratory load perceived by humans can be indirectly predicted by evaluating the variations in the nasal resistance and wall shear stress.28,58,59 As MRV measures only the velocity field and not pressure, the nasal resistance cannot be derived from the experimental data. In addition, calculating the wall shear stress near the wall is challenging owing to the partial volume effect and low resolution in the MRV data. Therefore, the validated CFD is used to obtain nasal resistance and wall shear stress. As discussed in Sec. III C, the laminar model is suitable under resting breathing conditions. To calculate the flow resistance of actual breathing, we performed CFD with air as the working fluid using the same grid. It should be noted that although air is a gaseous substance, it is considered an incompressible fluid at these low velocity conditions, and therefore, the overall fluid dynamics are similar to water which is an incompressible liquid, when the Reynolds number is matched for the two cases.

The overall nasal resistance can be expressed as Roverall = ΔPoverall/Qtotal, where ΔPoverall denotes the overall pressure drop from the nostril to nasopharynx and Qtotal indicates the respiratory flow rate. The overall pressure drop is depicted in Fig. 16, along with the three separate local pressure drops in the nasal valve, turbinate, and nasopharynx region, with volume rendering of equally spaced transparent slices of the overlapping pressure field. The pre-operative inhalation, post-operative inhalation, pre-operative exhalation, and post-operative exhalation cases are portrayed in Figs. 16(a)–16(d), respectively. In resting conditions with a flow rate of 100 ml/s, the overall nasal resistance before and after surgery was 0.047 and 0.016 Pa·s/ml for the inhalation peak, respectively, and 0.035 and 0.015 Pa·s/ml for the exhalation peak, respectively.60 After operation, the nasal resistance was reduced by 66% and 57% for the inhalation and exhalation peaks, respectively.

FIG. 16.

Local and overall pressure drop in the nasal valve, turbinate, and nasopharynx regions. (a) Pre-operative inhalation. (b) Post-operative inhalation. (c) Pre-operative exhalation. (d) Post-operative exhalation.

FIG. 16.

Local and overall pressure drop in the nasal valve, turbinate, and nasopharynx regions. (a) Pre-operative inhalation. (b) Post-operative inhalation. (c) Pre-operative exhalation. (d) Post-operative exhalation.

Close modal

In prior research, the nasal valves24,28,45,53 and nasopharynx52,53 were identified as regions of high local resistance, due to rapid cross-sectional contraction. However, the largest pressure drop actually occurs in the turbinate region. Since it occurs over a significantly longer length compared to the nasal valve and nasopharynx, the turbinate has received relatively less interest. The main topic of this study, though, is the improvement of flow characteristics before and after surgery for deviated nasal septum (DNS) or chronic hypertrophic rhinitis (CHR). Since these surgical procedures affect the turbinate region, the analysis of pressure drop in this area is necessary.

The nasal cavity is first segmented into three areas: the nasal vestibule region (from nostril to nasal valve), turbinate region (from nasal valve to choana), and nasopharynx (from choana to pharynx). The ratio of local to overall resistance is calculated and provided in Table III. In the pre-operative case, the highest resistance was observed in the turbinate region instead of the nasal valve region. As the patient suffered from DNS and CHR, the increase in the cross-sectional area was relatively small in the turbinate region, representing the air passage after the nasal valve, and the cross-sectional area/perimeter was significantly reduced (Fig. 17). This implied that the airflow passing through the long and narrow passage was affected by the wall in the turbinate region. In contrast, the cross-sectional area of the turbinate region increased remarkably after surgery, resulting in a cross-sectional area/perimeter increase by more than 45% compared to the pre-operation case. Consequently, the proportion of Rturbinate decreased by 26% for inhalation and 23% for exhalation after septoplasty.

TABLE III.

Fraction of local to overall nasal resistance.

F nasal valve F turbinate F nasopharynx F overall
Inhalation  Pre-operative  0.33  0.55  0.11 
Post-operative  0.42  0.29  0.29 
Exhalation  Pre-operative  0.21  0.71  0.08 
Post-operative  0.34  0.48  0.18 
F nasal valve F turbinate F nasopharynx F overall
Inhalation  Pre-operative  0.33  0.55  0.11 
Post-operative  0.42  0.29  0.29 
Exhalation  Pre-operative  0.21  0.71  0.08 
Post-operative  0.34  0.48  0.18 
FIG. 17.

Variations in the cross-sectional area along the air passage. (a) Pressure measurement plane in the air passage. (b) Local cross-sectional area and cross-sectional area/perimeter of pre- and post-operation cases.

FIG. 17.

Variations in the cross-sectional area along the air passage. (a) Pressure measurement plane in the air passage. (b) Local cross-sectional area and cross-sectional area/perimeter of pre- and post-operation cases.

Close modal
The left–right difference in nasal resistance can be obtained using the circuit model as follows in Strien et al. (Fig. 18).53  Q and R denote the flow rate and nasal resistance described above, and the subscripts denote the position and the left and right sides of the nasal cavity,
R overall = R left R right R left + R right + R pharynx ,
(5)
R left = R nasal valve , left + R turbinate , left ,
(6)
R right = R nasal valve , right + R turbinate , right ,
(7)
R left = Δ P nasal valve , left + Δ P turbinate , left Q left ,
(8)
R right = Δ P nasal valve , right + Δ P turbinate , right Q right .
(9)
FIG. 18.

Circuit model of nasal resistance during inhalation.

FIG. 18.

Circuit model of nasal resistance during inhalation.

Close modal

In the case of inhalation, the inlet pressure in the nostril is the same on the left and right sides (i.e., atmospheric pressure), and it also becomes the same after the two passages merge at the choana. Thus, the overall pressure loss between the nostril and choana is the same for both passages, and the nasal resistance becomes different between the left and right passages due to its respective flow rate. In the case of exhalation, the left and right exit pressures of the nostril remain the same (at atmospheric pressure), but the direction of flow is opposite to that of inhalation. The left and right flow rates and respective nasal resistance are listed in Table IV. For both inhalation and exhalation before surgery, the left/right flow rate ratio was approximately 1/2, resulting in a resistance ratio of approximately 2. However, after septoplasty, the difference between left and right became nearly negligible for both the flow rate and resistance, and the absolute value of the total resistance also decreased significantly.

TABLE IV.

Left and right respiratory flow rate and nasal resistance before and after operation.

Flow rate (ml/s) Nasal resistance (Pa·s/ml)
Qleft Qright Qleft/Qright Rleft Rright Rleft/Rright
Inhalation  Pre-operative  32.2  67.8  0.48  0.130  0.062  2.10 
Post-operative  51.5  48.5  1.06  0.022  0.023  0.94 
Exhalation  Pre-operative  31.3  68.7  0.45  0.144  0.065  2.20 
Post-operative  49.1  50.9  0.96  0.025  0.024  1.04 
Flow rate (ml/s) Nasal resistance (Pa·s/ml)
Qleft Qright Qleft/Qright Rleft Rright Rleft/Rright
Inhalation  Pre-operative  32.2  67.8  0.48  0.130  0.062  2.10 
Post-operative  51.5  48.5  1.06  0.022  0.023  0.94 
Exhalation  Pre-operative  31.3  68.7  0.45  0.144  0.065  2.20 
Post-operative  49.1  50.9  0.96  0.025  0.024  1.04 

The wall shear stress represents the velocity gradient near the wall, defined as τwss = μ · ∂u*/∂y*, where μ denotes the dynamic viscosity, u* indicates the velocity parallel to the wall, and y* denotes the distance in the normal direction of the wall. The calculated wall shear stress range was confirmed to have a similar range to that of previous studies.28 Comparing wall shear stress for the pre- and post-operative cases, strong wall shear is observed in the nasal valve and the lower portion of the middle turbinate in the pre-operative cases, as illustrated in Figs. 19(a) and 19(c). As previously shown in Fig. 17(b), when airflow passes from the nostril to the nasal valve, the area/perimeter decreases more significantly in the pre-operative case than the post-operative case. Due to the fast airflow through the nasal valve and the no-slip condition at the wall, the velocity gradient at the wall is large, resulting in a large shear stress. In the pre-operative case, there is a large deviation of nasal resistance between the left and right sides. A large amount of flow passes through the right side with less nasal resistance. The air passing through the narrow passage of the nasal valve enters the middle turbinate with significant inertia. The velocity gradient in the middle turbinate is also large due to the rapid airflow through the narrow passage, resulting in a large wall shear stress. However, as the narrow nasal cavity between the turbinates was widened in the post-operative case, the wall shear stress was significantly reduced compared to the pre-operative case, as depicted in Figs. 19(b) and 19(d).

FIG. 19.

Wall shear stress of nasal valve and lower part of middle turbinate (note the different orientation in this figure to show both the left nasal valve and the right middle turbinate). (a) Pre-operative inhalation. (b) Post-operative inhalation. (c) Pre-operative exhalation. (d) Post-operative exhalation.

FIG. 19.

Wall shear stress of nasal valve and lower part of middle turbinate (note the different orientation in this figure to show both the left nasal valve and the right middle turbinate). (a) Pre-operative inhalation. (b) Post-operative inhalation. (c) Pre-operative exhalation. (d) Post-operative exhalation.

Close modal

Although various functions of the nasal cavity are related to airflow, its noninvasive measurement in an experiment is difficult owing to the narrow and complex geometry of the nasal cavity. Previous noninvasive measurement attempts using 2D PIV could only examine a limited region within the nasal cavity at a fixed point in time. Therefore, this study employed noninvasive time-resolved MRV, which enables 3D three-component velocity measurements over the entire nasal cavity for each phase in the respiratory cycle. Using this technique, the total flow volume, cross-sectional area, and velocity field in the entire nasal cavity were measured at resting breathing condition for a patient suffering from a deviated nasal septum (DNS) and chronic hypertrophic rhinitis (CHR). The MRV data confirmed the increase in the nasal cross-sectional area after turbinoplasty, resulting in reduced nasal resistance and wall shear stress. A uniform left/right nasal area ratio was also achieved with septoplasty, resulting in a uniform left/right flow rate ratio wherein the deviation from unity was less than 0.8%. Consequently, the patient experienced more comfortable breathing after turbinoplasty and septoplasty.

The MRV data were directly compared with RANS CFD data to determine the most suitable CFD model (i.e., laminar, k–ε, k–ω, and SST) for nasal cavity flow. Similar to previous CFD studies of the nasal cavity, a steady-state flow within the respiration cycle was used.39 The E parameter from CFD was quite similar to that of the experimental results throughout the nasal cavity, suggesting that this method of comparing transient experimental data with steady CFD data is valid. Despite the low Reynolds number of 591 at the peak respiration point in the trachea, we employed different RANS models considering the possible occurrence of local turbulence owing to the complex geometry and rapid variations in the cross-sectional area. The results confirmed that the laminar, k–ω, and SST models were practically identical to each other and also quite similar to the experiment in all regions of the nasal cavity. This implies that the Reynolds number is low enough such that the flow can be considered fully laminar. Therefore, considering the cost of simulation, the laminar model can be reasonably applied for resting breathing conditions (100 ml/s). However, further study is required to determine the variation across different patients or to experiment on models representative of human nasal geometry.41,46

In the future, we intend to conduct similar studies under conditions of increased intensity breathing, such as moderate and extreme exercise, sniffing, coughing, and sneezing, using 4D MRV. We also plan to validate additional RANS models such as the low Reynolds number (LRN) k–ω and transition SST models, while also comparing results to large eddy simulation (LES). Conducting experiments at steady flow conditions will also allow us to examine the effect of higher Reynolds numbers where turbulence might play a stronger role.

This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (Ministry of Education) (No. NRF-2021R1A6A3A13040106) and also supported by the Interdisciplinary Research Initiatives Program from College of Engineering and College of Medicine, Institute of Engineering Research, Seoul National University (Grant No. 800-20200285, 2020-IOER-0), SNUH Research Fund (Grant No. 0320220320), and Institute of Advanced Machines and Design at Seoul National University.

The authors have no conflicts to disclose.

Kyuho Han and Sung-Gwang Lee contributed equally to this work.

Kyuho Han: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Funding acquisition (equal); Investigation (equal); Methodology (equal); Project administration (equal); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Sung-Gwang Lee: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (equal); Methodology (equal); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Kwanwoo Kim: Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Validation (supporting); Visualization (supporting); Writing – original draft (supporting). Baren Jeong: Data curation (supporting); Formal analysis (supporting); Funding acquisition (supporting); Methodology (supporting); Project administration (equal); Resources (supporting); Visualization (supporting); Writing – original draft (supporting). Munyoung Paek: Methodology (supporting); Software (supporting); Writing – original draft (supporting). Whal Lee: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Resources (equal); Supervision (equal); Writing – original draft (supporting). Wontae Hwang: Conceptualization (equal); Funding acquisition (equal); Project administration (equal); Resources (equal); Supervision (lead); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.

1.
S. E.
Churchill
,
L. L.
Shackelford
,
J. N.
Georgi
, and
M. T.
Black
, “
Morphological variation and airflow dynamics in the human nose
,”
Am. J. Hum. Biol.
16
(
6
),
625
638
(
2004
).
2.
A.
Thomson
and
L. D.
Buxton
, “
Man’s nasal index in relation to certain climatic conditions
,”
J. R. Anthropol. Inst. G. B. Irel.
53
,
92
122
(
1923
).
3.
M. H.
Wolpoff
, “
Climatic influence on the skeletal nasal aperture
,”
Am. J. Phys. Anthropol.
29
(
3
),
405
423
(
1968
).
4.
R. G.
Franciscus
and
J. C.
Long
, “
Variation in human nasal height and breadth
,”
Am. J. Phys. Anthropol.
85
(
4
),
419
427
(
1991
).
5.
K.
Liener
,
R.
Leiacker
,
J.
Lidemann
,
G.
Rettinger
, and
T.
Keck
, “
Nasal mucosal temperature after exposure to cold, dry air and hot, humid air
,”
Acta Oto-Laryngol.
123
(
7
),
851
856
(
2003
).
6.
D. F.
Proctor
and
I. H. P.
Andersen
,
The Nose: Upper Airway Physiology and the Atmospheric Environment
(
Elsevier Biomedical Press
,
Amsterdam
,
1982
), pp.
23
43
.
7.
C. S.
Kim
,
B. K.
Moon
,
D. H.
Jung
, and
Y. G.
Min
, “
Correlation between nasal obstruction symptoms and objective parameters of acoustic rhinometry and rhinomanometry
,”
Auris, Nasus, Larynx
25
(
1
),
45
48
(
1998
).
8.
K.
Wiesmiller
,
T.
Keck
,
G.
Rettinger
,
R.
Leiacker
,
R.
Dzida
, and
J.
Lindemann
, “
Nasal air conditioning in patients before and after septoplasty with bilateral turbinoplasty
,”
Laryngoscope
116
(
6
),
890
894
(
2006
).
9.
A. W.
Proetz
, “
XLI air currents in the upper respiratory tract and their clinical importance
,”
Ann. Otol., Rhinol. Laryngol.
60
(
2
),
439
467
(
1951
).
10.
M.
Stuiver
,
Biophysics of the Sense of Smell
(
Uitgeverij Excelsior
,
1958
).
11.
H.
Masing
, “
Investigations about the course of flow in the nose model
,”
Arch. Klin. Exp. Ohr-, Nas- U Kehlk. Heilk.
189
,
371
381
(
1967
).
12.
D. J.
Brain
,
D. F.
Proctor
, and
L. M.
Reid
,
Respiratory Defense Mechanisms, Part I
(
Dekker
,
New York
,
1977
).
13.
M.
Girardin
,
E.
Bilgen
, and
P.
Arbour
, “
Experimental study of velocity fields in a human nasal fossa by laser anemometry
,”
Ann. Otol. Rhinol. Laryngol.
92
(
3
),
231
236
(
1983
).
14.
S.
Schreck
,
K. J.
Sullivan
,
C. M.
Ho
, and
H. K.
Chang
, “
Correlations between flow resistance and geometry in a model of the human nose
,”
J. Appl. Physiol.
75
(
4
),
1767
1775
(
1993
).
15.
I.
Hahn
,
P. W.
Scherer
, and
M. M.
Mozell
, “
Velocity profiles measured for airflow through a large-scale model of the human nasal cavity
,”
J. Appl. Physiol.
75
(
5
),
2273
2287
(
1993
).
16.
L. M.
Hopkins
,
J. T.
Kelly
,
A. S.
Wexler
, and
A. K.
Prasad
, “
Particle image velocimetry measurements in complex geometries
,”
Exp. Fluids
29
(
1
),
91
95
(
2000
).
17.
J. T.
Kelly
,
A. K.
Prasad
, and
A. S.
Wexler
, “
Detailed flow patterns in the nasal cavity
,”
J. Appl. Physiol.
89
(
1
),
323
337
(
2000
).
18.
S. K.
Kim
and
S. K.
Chung
, “
An investigation on airflow in disordered nasal cavity and its corrected models by tomographic PIV
,”
Meas. Sci. Technol.
15
(
6
),
1090
(
2004
).
19.
S. K.
Chung
,
Y. R.
Son
,
S. J.
Shin
, and
S. K.
Kim
, “
Nasal airflow during respiratory cycle
,”
Am. J. Rhinol.
20
(
4
),
379
384
(
2006
).
20.
T.
Soodt
,
F.
Schröder
,
M.
Klaas
,
T.
van Overbrüggen
, and
W.
Schröder
, “
Experimental investigation of the transitional bronchial velocity distribution using stereo scanning PIV
,”
Exp. Fluids
52
(
3
),
709
718
(
2012
).
21.
E.
Nof
,
S.
Bhardwaj
,
P.
Koullapis
,
R.
Bessler
,
S.
Kassinos
, and
J.
Sznitman
, “
In vitro–in silico correlation of three-dimensional turbulent flows in an idealized mouth-throat model
,”
PLoS Comput. Biol.
19
(
3
),
e1010537
(
2023
).
22.
S. K.
Kim
,
Y.
Na
,
J. I.
Kim
, and
S. K.
Chung
, “
Patient specific CFD models of nasal airflow: Overview of methods and challenges
,”
J. Biomech.
46
(
2
),
299
306
(
2013
).
23.
J.
Lindemann
,
H. J.
Brambs
,
T.
Keck
,
K. M.
Wiesmiller
,
G.
Rettinger
, and
D.
Pless
, “
Numerical simulation of intranasal airflow after radical sinus surgery
,”
Am. J. Otolaryngol.
26
(
3
),
175
180
(
2005
).
24.
S.
Ozlugedik
,
G.
Nakiboglu
,
C.
Sert
,
A.
Elhan
,
E.
Tonuk
,
S.
Akyar
, and
I.
Tekdemir
, “
Numerical study of the aerodynamic effects of septoplasty and partial lateral turbinectomy
,”
Laryngoscope
118
(
2
),
330
334
(
2008
).
25.
X. B.
Chen
,
H. P.
Lee
,
V. F. H.
Chong
, and
D. Y.
Wang
, “
A computational fluid dynamics model for drug delivery in a nasal cavity with inferior turbinate hypertrophy
,”
J. Aerosol Med. Pulm. Drug Delivery
23
(
5
),
329
338
(
2010
).
26.
G. J.
Garcia
,
J. S.
Rhee
,
B. A.
Senior
, and
J. S.
Kimbell
, “
Septal deviation and nasal resistance: An investigation using virtual surgery and computational fluid dynamics
,”
Am. J. Rhinol. Allergy
24
(
1
),
e46
e53
(
2010
).
27.
J. S.
Rhee
,
S. S.
Pawar
,
G. J.
Garcia
, and
J. S.
Kimbell
, “
Toward personalized nasal surgery using computational fluid dynamics
,”
Arch. Facial Plast. Surg.
13
(
5
),
305
310
(
2011
).
28.
Y.
Na
,
K. S.
Chung
,
S. K.
Chung
, and
S. K.
Kim
, “
Effects of single-sided inferior turbinectomy on nasal function and airflow characteristics
,”
Respir. Physiol. Neurobiol.
180
(
2–3
),
289
297
(
2012
).
29.
Z.
Zhang
and
C.
Kleinstreuer
, “
Low-Reynolds-number turbulent flows in locally constricted conduits: A comparison study
,”
AIAA J.
41
(
5
),
831
840
(
2003
).
30.
K.
Bass
and
P. W.
Longest
, “
Recommendations for simulating microparticle deposition at conditions similar to the upper airways with two-equation turbulence models
,”
J. Aerosol Sci.
119
,
31
50
(
2018
).
31.
X.
Chen
,
Y.
Feng
,
W.
Zhong
, and
C.
Kleinstreuer
, “
Numerical investigation of the interaction, transport and deposition of multicomponent droplets in a simple mouth-throat model
,”
J. Aerosol Sci.
105
,
108
127
(
2017
).
32.
X.
Chen
,
X.
Zhou
,
X.
Xia
,
X.
Xie
,
P.
Lu
, and
Y.
Feng
, “
Modeling of the transport, hygroscopic growth, and deposition of multi-component droplets in a simplified airway with realistic thermal boundary conditions
,”
J. Aerosol Sci.
151
,
105626
(
2021
).
33.
C. J.
Elkins
and
M. T.
Alley
, “
Magnetic resonance velocimetry: Applications of magnetic resonance imaging in the measurement of fluid motion
,”
Exp. Fluids
43
(
6
),
823
858
(
2007
).
34.
S.
Baek
,
S.
Lee
,
W.
Hwang
, and
J. S.
Park
, “
Experimental and numerical investigation of the flow in a trailing edge ribbed internal cooling passage
,”
J. Turbomach.
141
(
1
),
011012
(
2019
).
35.
C. J.
Elkins
,
M.
Markl
,
N.
Pelc
, and
J. K.
Eaton
, “
4D Magnetic resonance velocimetry for mean velocity measurements in complex turbulent flows
,”
Exp. Fluids
34
(
4
),
494
503
(
2003
).
36.
M.
Markl
,
A.
Harloff
,
T. A.
Bley
,
M.
Zaitsev
,
B.
Jung
,
E.
Weigang
,
M.
Langer
,
J.
Hennig
, and
A.
Frydrychowicz
, “
Time‐resolved 3D MR velocity mapping at 3T: Improved navigator‐gated assessment of vascular anatomy and blood flow
,”
Magn. Reson. Imaging
25
(
4
),
824
831
(
2007
).
37.
A. J.
Banko
,
F.
Coletti
,
C. J.
Elkins
, and
J. K.
Eaton
, “
Oscillatory flow in the human airways from the mouth through several bronchial generations
,”
Int. J. Heat Fluid Flow
61
,
45
57
(
2016
).
38.
R.
Medero
,
C.
Hoffman
, and
A.
Roldán-Alzate
, “
Comparison of 4D flow MRI and particle image velocimetry using an in vitro carotid bifurcation model
,”
Ann. Biomed. Eng.
46
(
12
),
2112
2122
(
2018
).
39.
S.
Jalal
,
T.
van de Moortele
,
O.
Amili
, and
F.
Coletti
, “
Steady and oscillatory flow in the human bronchial tree
,”
Phys. Rev. Fluids
5
(
6
),
063101
(
2020
).
40.
C. J.
Collier
,
L.
Langlois
,
Y.
Ow
,
C.
Johansson
,
M.
Giammusso
,
M. P.
Adams
,
K. R.
O'Brien
, and
S.
Uthicke
, “
Losing a winner: Thermal stress and local pressures outweigh the positive effects of ocean acidification for tropical seagrasses
,”
New Phytol.
219
(
3
),
1005
1017
(
2018
).
41.
D. D.
Borup
,
L. E.
Engel
,
C. J.
Elkins
, and
J. K.
Eaton
, “
Transport and dispersion of particle-Laden streaks in a standardized human nasal geometry
,”
Exp. Fluids
61
(
2
),
43
(
2020
).
42.
C. J. T.
Spence
,
N. A.
Buchmann
,
M. C.
Jermy
, and
S. M.
Moore
, “
Stereoscopic PIV measurements of flow in the nasal cavity with high flow therapy
,”
Exp. Fluids
50
(
4
),
1005
1017
(
2011
).
43.
J. H.
Zhu
,
H. P.
Lee
,
K. M.
Lim
, and
S. J.
Lee
, “
Evaluation and comparison of nasal airway flow patterns among three subjects from Caucasian, Chinese and Indian ethnic groups using computational fluid dynamics simulation
,”
Respir. Physiol. Neurobiol.
175
(
1
),
62
69
(
2011
).
44.
X. B.
Chen
,
S. C.
Leong
,
H. P.
Lee
,
V. F.
Chong
, and
D. Y.
Wang
, “
Aerodynamic effects of inferior turbinate surgery on nasal airflow—a computational fluid dynamics model
,”
Rhinology
48
(
4
),
394
400
(
2010
).
45.
J.
Wen
,
K.
Inthavong
,
J.
Tu
, and
S.
Wang
, “
Numerical simulations for detailed airflow dynamics in a human nasal cavity
,”
Respir. Physiol. Neurobiol.
161
(
2
),
125
135
(
2008
).
46.
J.
Choi
,
M. H.
Tawhai
,
E. A.
Hoffman
, and
C. L.
Lin
, “
On intra- and intersubject variabilities of airflow in the human lungs
,”
Phys. Fluids
21
(
10
),
101901
(
2009
).
47.
D. L.
Jan
,
A. H.
Shapiro
, and
R. D.
Kamm
, “
Some features of oscillatory flow in a model bifurcation
,”
J. Appl. Physiol.
67
(
1
),
147
159
(
1989
).
48.
M.
Bruschewski
,
D.
Freudenhammer
,
W. B.
Buchenberg
,
H. P.
Schiffer
, and
S.
Grundmann
, “
Estimation of the measurement uncertainty in magnetic resonance velocimetry based on statistical models
,”
Exp. Fluids
57
(
5
),
83
(
2016
).
49.
T. L.
Chan
,
R. M.
Schreck
, and
M.
Lippmann
, “
Effect of the laryngeal jet on particle deposition in the human trachea and upper bronchial airways
,”
J. Aerosol Sci.
11
(
5–6
),
447
459
(
1980
).
50.
G.
Mylavarapu
,
S.
Murugappan
,
M.
Mihaescu
,
M.
Kalra
,
S.
Khosla
, and
E.
Gutmark
, “
Validation of computational fluid dynamics methodology used for human upper airway flow simulations
,”
J. Biomech.
42
(
10
),
1553
1559
(
2009
).
51.
H.
Shi
,
C.
Kleinstreuer
, and
Z.
Zhang
, “
Dilute suspension flow with nanoparticle deposition in a representative nasal airway model
,”
Phys. Fluids
20
(
1
),
013301
(
2008
).
52.
C.
Kleinstreuer
and
Z.
Zhang
, “
Airflow and particle transport in the human respiratory system
,”
Annu. Rev. Fluid Mech.
42
,
301
334
(
2010
).
53.
J.
Van Strien
,
K.
Shrestha
,
S.
Gabriel
,
P.
Lappas
,
D. F.
Fletcher
,
N.
Singh
, and
K.
Inthavong
, “
Pressure distribution and flow dynamics in a nasal airway using a scale resolving simulation
,”
Phys. Fluids
33
(
1
),
011907
(
2021
).
54.
C.
Li
,
J.
Jiang
,
H.
Dong
, and
K.
Zhao
, “
Computational modeling and validation of human nasal airflow under various breathing conditions
,”
J. Biomech.
64
,
59
68
(
2017
).
55.
D. J.
Doorly
,
D. J.
Taylor
, and
R. C.
Schroter
, “
Mechanics of airflow in the human nasal airways
,”
Respir. Physiol. Neurobiol.
163
(
1–3
),
100
110
(
2008
).
56.
D.
Doorly
,
D. J.
Taylor
,
P.
Franke
, and
R. C.
Schroter
, “
Experimental investigation of nasal airflow
,”
Proc. Inst. Mech. Eng., Part H
222
(
4
),
439
453
(
2008
).
57.
A. M.
Padilla
,
The Effect of Upstream Perturbations on 3D Annular Diffusers
(
Stanford University
,
2012
).
58.
A. H.
Suzina
,
M.
Hamzah
, and
A. R.
Samsudin
, “
Objective assessment of nasal resistance in patients with nasal disease
,”
J. Laryngol. Otol.
117
(
8
),
609
613
(
2003
).
59.
D.
Elad
,
S.
Naftali
,
M.
Rosenfeld
, and
M.
Wolf
, “
Physical stresses at the air-wall interface of the human nasal cavity during breathing
,”
J. Appl. Physiol.
100
(
3
),
1003
1010
(
2006
).
60.
A. A.
Borojeni
,
G. J.
Garcia
,
M. G.
Moghaddam
,
D. O.
Frank-Ito
,
J. S.
Kimbell
,
P. W.
Laud
,
L. J.
Koenig
, and
J. S.
Rhee
, “
Normative ranges of nasal airflow variables in healthy adults
,”
Int. J. Comput. Assisted Radiol. Surg.
15
(
1
),
87
98
(
2020
).