Investigation of the effects of nasal surgery on nasal cavity flow using magnetic resonance velocimetry and computational fluid dynamics

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.¹ 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.⁵ The airflow also becomes more turbulent in the nostril as it is directed downward through a tortuous path.⁶ 

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.⁷ Wiesmiller et al. investigated the effect of surgical improvement based on in vivo measurements of nasal temperature and humidity before and after septoplasty.⁸ 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.¹² Thereafter, Girardin et al. measured the quantitative velocity field using laser doppler velocimetry.¹³ Schreck et al. conducted pressure measurements, hot-wire anemometry, and flow visualization using dye in water.¹⁴ Hahn et al. fabricated a right nasal model of 1:20 enlarged scale to measure the velocity profile using a hot-film anemometer probe.¹⁵ 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.²² Several CFD studies discussed the correlation between flow and pathology in post-operation nasal models.^23–28 Ozlugedik et al.²⁴ 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.²⁸ 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 Kleinstreuer²⁹ 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 Lingest³⁰ 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.³⁵ Therefore, this technique has been applied to biofluid measurements within vascular and respiratory systems.^36–38 Banko et al.³⁷ examined the flow through the larynx, trachea, and bronchus under oral breathing conditions. Jalal et al.³⁹ used 4D MRV to investigate the velocity field from the middle of the trachea to the bronchus. Collier et al.⁴⁰ utilized MRV as a validation tool for CFD simulations of the 3D mean flow in the central airway tree. Borup et al.⁴¹ 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.

The reconstructed full-scale test section is illustrated in Fig. 2. The nasal cavity model was constructed using 3D Builder.⁴² 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.⁴³ 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.

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 mm³ in the x-, y-, and z-directions, using a voxel size of 0.5 × 0.5 × 1.0 mm³. 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.

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.

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 × 10⁶ meshes in the pre-operative model and 11.5 × 10⁶ 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.

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.⁵⁰ 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.

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.

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.

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.⁵¹ 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.

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.

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.⁵² 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.²⁸ 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.⁵⁴ 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).

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.

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 R_overall = ΔP_overall/Q_total, where ΔP_overall denotes the overall pressure drop from the nostril to nasopharynx and Q_total 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.⁶⁰ After operation, the nasal resistance was reduced by 66% and 57% for the inhalation and exhalation peaks, respectively.

In prior research, the nasal valves^24,28,45,53 and nasopharynx^52,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 R_turbinate decreased by 26% for inhalation and 23% for exhalation after septoplasty.

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.

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.²⁸ 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).

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.³⁹ 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.