Background: The fluid flow in channels with deviation change in the cross-sectional area in a part of the channel can have a turbulent behavior. Nasal Septal Deviation (NSD) is such channel and can disturb airflow. Computational Fluid Dynamics (CFD) can show the critical points that cause turbulence or blockage of the flow by analyzing the airflow passing through the nasal cavity in detail.
Objective: We evaluated a patient with nasal septal deviation before and after correcting septum surgery (septoplasty). With the help of Computational Fluid Dynamics (CFD) by showing the patterns of airflow through the patient's nasal cavity and accurately locating the deviation of the septum along with showing the disturbance of the passing airflow to the surgeon, it is possible to preserve healthy tissues. 
Methods: First, Computed Tomography (CT) of the nasal cavity and paranasal sinuses was done. Then, using Mimics software, the basic 3D geometry of the nose was made. Due to the point-like nature of the created shape, we converted the 3D geometry from the point to a linear geometry and finally to a volumetric geometry by CATIA software. Considering the transient flow from the nasal tip, simulation of the numerical method at air flow rates of 15, 17.4, 20, 25 and 40 lit/min was done using the k-ω SST turbulence model. We considered the non-slip and smooth nasal wall and steady flow at 22 ᵒC.  According to CFD data and turbulence area the exact location of nasal cavity which caused obstruction was shown to surgeon and septoplasty was done. 3D modeling using CFD was done and airflow parameters were evaluated again three months after surgery to evaluate development of airflow behavior.
Results: The peak velocity was found at 3 to 6 cm from the nostril before septoplasty, but the peak velocity was at the nasal valve after septoplasty. The airflow resistance decreased from 0.97 Pa×min/lit to 0.53 Pa×min/lit after septoplasty for a flow rate of 17.4 lit/min (the beginning of airflow disturbance). Also, the pressure drop for a flow rate of 17.4 liters per minute was 16.878 Pa, which decreased to 9.222 Pa after surgery. The highest wall shear stress was also at a distance of 4.1 and 5.07cm from nasal tip before and after septoplasty respectively. The highest wall shear stress was 3.85×10-2 Pa before septoplasty, but after septoplasty it decreased to 2.78×10-2 Pa.
Conclusion: Nasal septal deviation increases the velocity, resistance, pressure drop and creates a vortex in the airflow, which leads to obstruction of the airway passing through the nasal cavity. CFD showed difference in nasal airflow parameters in patients with deviated nasal septum before and after septoplasty, so it can help surgeon find the exact area which causes obstruction and septoplasty can be done according to CFD data.

The nasal cavity is known to be one of the most important parts of the human airway system, which provides the first protection stage of the lungs by heating and moistening the air [1–3]. The upper respiratory tract from the valve of the nose to the entrance of the larynx plays an important role in delivering oxygen to the lungs, leads to absorbing in blood. However, the cross-sectional regularity of the air passage from the nasal cavity to the larynx plays an important role in getting enough oxygen to the lungs and improving respiration. Hence, a better understanding of how streamlines of airflow, the velocity and pressure at obstructive and critical points in the nasal cavity on the inhaled airflow is essential [4, 5].
The complicated structure associated with the nasal anatomy makes it difficult to measure the resistance of the nose [6]. Determining the nasal anatomy and physiology is difficult due to the small size of these airways. Several researchers have performed studies of airflow in the nasal cavity using measuring devices such rhinomanometry and acoustic rhinometry [7–10].
Harkma et al. and Grotberg reported that determining the exact location of the obstruction in the nasal cavity is the first step in identifying the most important factors playing a role as understanding nasal and upper airway pathology [11]. Airflow through the human upper respiratory tract has been studied numerically and experimentally by several researchers [12–26].
Rhinomanometry is used to measure the pressure of airflow through the nasal airway and acoustic rhinometry is used to measure the cross section of the airway in different planes of the nose. However, measuring the exact velocity of the airflow and assessing the cross-sectional area in each part of the nasal cavity is difficult [27–31].
The reliability of the results obtained using this device depends on the favorable cooperation of the operator of the rhinomanometric device and the patient, the calibration of the device and correct executive instructions by the physician and the use of standard techniques. Due to the inherent limitations of existing measuring instruments, Computational Fluid Dynamics (CFD) has been proposed as a suitable option [32].
CFD which refers to the use of numerical methods to solve the Partial Differential Equations (PDE) governing fluid flow, is becoming an increasingly popular research tool in the study and prediction of fluid dynamics [33]. 
Tretiakow, D et al. used nasal CT scans from 16 patients to process images and made a model of the human nasal cavity and sinuses using a 3D slicer for segmentation and new surface model generation. Further processing was performed in Autodesk® Mesh Mixer TM (https://meshmixer.com/). The underlying equations were discretized using the finite volume method. The corresponding system of algebraic equations was then solved using Open foam software (https://www.openfoam.com/). The CFD results were presented using sample 3D models of patient 1 (normal) and patient 2 (pathological changes). As a conclusion, CFD is recommended to be used in clinical practice and also in research as a beneficial objective tool to evaluate the nasal airway [34].
CFD is also useful in evaluating the outcomes of nasal surgery. For example, partial or complete inferior turbinate reduction is assessed and compared in Lee et al.’s study using CFD analysis [35]. 
The mechanism of feeling obstruction in a patient with atrophic rhinitis (excessive dilation of nasal cavity) using CFD is discussed well in Ya Zhang et al.’s study [36]. 
In this study, we aimed at using CFD data for determining the exact area that caused obstruction and suggesting the surgeon to do septoplasty according to mentioned area not only saves the healthy nasal tissue, but also improves the results of septoplasty. 

Reconstruction of geometry

The patient was a 45-year-old man of Iranian nationality who had breathing problems due to obstruction of the right nasal cavity. The NOSE questionnaire score for the patient was 12 (the NOSE questionnaire is uploaded as a supplementary material). The patient has not had a common cold, sinusitis, nasal polyps, previous surgery on the nose or trauma to the face in the past six months. First, Computed Tomography (CT scan) images of the face were performed by a radiologist at Ghaem Hospital, Mashhad University of Medical Sciences, Iran (figure 1).

Figure 1: Tomography images taken from the patient before septoplasty

In this research, we used version 19 of Mimics (Materialize Mimics) software. CT scan images were imported into Mimics medical imaging software (Materialsize, Plymouth, MI, USA). The 3D reconstruction of the 

nasal cavity was performed by removing the paranasal sinuses due to the prevention of tube airflow disturbance in the tube-like channel of the nasal cavity (Figure 2).

Figure 2: Reconstructed images of the patient's nasal cavity before surgery by Mimics software

Mesh geometry
The output of Mimics software is in the form of discontinuous lines named. STL format. This format cannot be used in Ansys Meshing and Fluent software, so the output file of Mimcs 3D images was transferred to Catia software. With the formation of cloud points in the Catia software, the 3D point image was converted into a linear image, and then the volume model was developed. After correcting the smoothness of the surfaces and removing the extra points, the boundary conditions were determined. Output file from the Catia software was generated as. STEP (Figure 3).

Figure 3: Modified images of nasal cavity in Catia software

Finally, the 3D geometry of the nasal cavity was imported into the Ansys ICEM CFD software as a .STEP file, and the meshing of the 3D geometry was generated. The initial mesh for the case before nasal cavity surgery was 1213339. By checking the meshes, we made the grid cells smaller and increased the number of meshes for the nasal cavity. By increasing the number of meshes from 4607378 cells to 5673085 cells, there is no significant change in the average pressure and average wall shear stress.
After our interpretation and numerical analysis with the help of CFD on the airflow passing through the patient's nasal cavity and determining the critical cross-sectional area that was the cause of flow disturbance and transient region. ENT surgeon done the septoplasty according to CFD data.
3 months after surgery, nasal cavity CT csan was repeated. Then the reconstruction steps of 3D images formation were done. Next, the airflow was analyzed using Ansys meshing and Fluent. At last, we investigated the change in nasal cavity airflow characteristics before and after septoplasty. Also, the NOSE questionnaire was used 3 months after septoplasty and it was 4.

Model of Air Flow in Nasal Cavity

The most successful way to discretize the flow problem is to use the control volume method, also known as the finite volume method or Euler method. As the name suggests, the problem is split into a control volume and the conservation laws of physics apply to this volume. The governing equation is:

Continuity equation (mass conservation);

Momentum equation (Newton's second law);

Reynolds’ transport theorem.

Reynolds' transport theorem is used to convert the governing equations to Euler form [37].

Conservation of mass means that the amount of mass remains constant over time. Therefore, if there is mass flow through a volume, the accumulation rate within that volume is equal to the net outflow. In vector notation:

(1)

This is also called the continuity equation for compressible fluids in unsteady flow. For an incompressible fluid, the continuity equation is:

(2)

Reynolds’ Transport Theory

The rate of change of the general property per unit volume can be expressed similarly to the rate of change of mass per unit volume (Equation (1)) as:

(3)

The characteristic rate of change is equal to the sum of the control volume rate of increase and the net outflow.

The Momentum Equation

Once Reynolds’ transport theory is established, it can be applied to Newton's second law to compute the rate of change of momentum using a control volume approach: 

(4)

Using:

	(5)
	(6)

The total force acting on a volume is the sum of pressure, gravity (and other body forces), shear force and normal stress [37]. Applying Reynolds' transport theorem and Newton's second law to a three-dimensional infinitesimal control volume yields the three-dimensional Navier–Stokes differential equation for laminar flow [38]:

)

(7)

Due to the turbulence in this case, velocity fluctuations and additional shear stresses, the so-called Reynolds stresses, must be taken into account. Velocity is modeled as the sum of mean and variation, u = U + u’. This produces a slightly different version of Navier–Stokes:

(8)

Using the Boussinesq approximation for the Reynolds stress term (turbulent stress) on the right-hand side produces:

(9)

After insertion and rearrangement, the Reynolds-averaged Navier–Stokes equations for turbulence are obtained [39].

For steady state, the transient term on the left is of course zero, so the convective term on the left is left alone. P is the pressure; k is the kinetic energy and the remaining two terms are diffusion (turbulent stress) and viscous stress. is the Kronecker delta required for the formula to give correct results for normal Reynolds stress. 

The turbulent kinetic energy equation is:

(10)

The pseudo-vortices equation () is:

(11)

Note that the addition notation i, j = 1, 2, 3 denotes the components of the velocity vector and spatial coordinates. Turbulent viscosities and their formulations are described in detail in ref. [39]. Note that the first terms in Equations (6)– (8) are ignored since steady-state conditions are applied in this study.

Numerical Solution

Combining the governing equations yields a closed system with an equal number of equations and unknowns. What now remains is a way to solve these equations numerically. The Ansys Fluent 18 solver for incompressible steady-state turbulence, called Simple Algorithm, uses the SIMPLE algorithm to solve the flow. 

The Simple Method

SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) starts with a (preferably qualified) guess of initial and boundary conditions for the pressure field. The guessed value is marked with an asterisk, . is used to the discretized momentum equations to calculate velocities, . These velocities are referred to and treated as guessed values.

The continuity equation is used to express a correction to the pressure, . The guessed and correction values are summed to get the actual pressure. Similarly, a velocity correction is added to :

	(12)
	(13)

After correcting for pressure and velocity, the other discretized transport equations are solved. This completes one iteration. The process is repeated using the corrected pressures and velocities as initial estimates for the next iteration. This continues until the convergence criterion is met. Different criteria may exist for different flow characteristics. The convergence criteria for this paper are described in Section 2.3.7. A flowchart of the SIMPLE method by Versteeg and Malalasekera is shown below [40].

There are several versions of this algorithm, for example, SIMPLE revised or the SIMPLER method, which directly computes the pressure field via the pressure discretization equation [40]. Nevertheless, only the original SIMPLE method was used in this study, as it forms the basis for the chosen software solver Simple.

Under-Relaxation

The SIMPLE method is prone to instability. One way to prevent this is by introducing under-relaxation coefficients,, for the pressure and velocity corrections. An under-relaxation coefficient is a factor between 0 and 1 that is used to dampen corrected values:

	(14)
	(15)

Even if the solution oscillates, the amplitude will be small, so the solution may not oscillate or diverge. However, this slows down convergence. The optimal relaxation factor depends on mesh, flow problem and iteration method and is difficult to predict [40].

Numerical Schemes

The most important basic properties of discretization methods are:

Conservatism: The same number of properties entering the cell area on one side exit the same area on the other side;

Boundedness: The property value of a point is within the range spanned by the boundary

Boundary   < value At Node < boundary

(16)

Transportiveness: The influence of the property ∅ on the value of a cell by neighboring cells depends on the convective–diffusion relationship in the flow. This relationship is called the Peclet number, Pe. The influence of neighboring cells should be influenced according to the Peclet number.

(17)

For pure diffusion, −Pe > 0; for pure convection, −Pe > 1;

To use routine comparison;

The purchase number scheme is similar to other purchases.

It is hard to complete all the processes using one product. The desired properties, such as second-order precision and boundedness, can be mutually exclusive. Quality also comes at a cost. Sometimes the best systems are too expensive. However, the scheme works well, even if it has only first-order precision or is less portable, as long as it is suitable for the problem at hand. For example, the central derivative scheme is sufficient for pure diffusion, but becomes less transportable as increases [40]. In this case, the scheme was chosen to achieve second-order accuracy. 

Boundary Conditions

To create a physical or computational model, you must identify limits and set appropriate values for those limits. The nasal airways have boundaries that change over time and are usually lined with a thin layer of mucus. The shape of the nose also affects the internal flow due to the inspiratory cycle and the flow entering the nostrils. During the respiratory cycle, the airway is imprinted under the nasopharynx. Due to the aforementioned complexity and the authors' knowledge, no model yet covers all these features. In the current study, three different constant inlet velocities of 15, 17.4 and 20 LPM were used as inlet conditions to simulate normal breathing in sleep, rest and relaxation situations. In addition, atmospheric pressure conditions are set at the outlet. As already mentioned, the effects of gravitational acceleration are ignored and a no-slip condition is used on the walls of the model. 

The initial values for k and ω at the inlet are determined using the following relations [41]:

	(18)
	(19)

where is the hydraulic diameter and I is the turbulence intensity obtained from the following equation [40]:

(20)

The nasal wall is considered rigid and non-slip [31,42]. A mass flow rate is applied to the nostril inlet, which corresponds to a specific air inlet (liters/minute). The flow range considered in this study is 15, 17.4 and 20 LPM and exists in the laminar and turbulent regimes for both cases [9,43,44]. As a general rule, for an adult nasal cavity, flow below 15 LPM is considered laminar, between 15 and 17.4 LPM is transient and above 17.4 LPM is turbulent. For the nasopharynx, the boundary condition “outflow” is considered, which assumes a fully developed flow. This simulation does not take into account the presence of mucus layers and nose hairs. The physical characteristics of the fluid (air) used in this calculation are a density of 1.225 and a kinematic viscosity of 1.7894 × .

The governing equations are discretized over the control volume. Integrating these discretized equations yields a set of equations in algebraic form. Simulations were performed using the CFD solver Ansys Fluent 18 R2. The SIMPLE algorithm was chosen to correlate velocity and pressure corrections. For higher accuracy, second-order schemes of momentum, turbulent kinetic energy and specific dissipation factor are used, and a convergence criterion of four orders of magnitude is adopted.

Solution Methods

In this study, Fluent 18 is applied and the convergence criterion is . The solution algorithms are shown in Figure 6.

The initial values ​​of k and ω in the input are determined using the following relationships [41].

	(18)
	(19)

Where is the hydraulic diameter and I is the turbulence intensity, which is obtained from the following equation.

(20)

The nasal wall is considered rigid and non-slip [30 and 42]. The mass flow rate applied to the nasal cavity inlet corresponds to a specific air inlet (liters per minute or LPM). The range of flow rate considered in this study is 15, 17.4 and 20 LPM and it exists in transient and turbulent regime for both cases [9-43-44]. As a general rule, for the nasal cavity of adults, a flow of less than 15 LPM is considered calm, between 15 and 17.4 LPM transient, and above 17.4 LPM turbulent. For the exit of the nasal cavity (nasopharynx), the "exit" boundary condition is considered, assuming a fully developed flow. In this simulation, we do not consider the existence of mucous layers and nasal villi. The physical characteristics of the fluid (air) used in this calculation are density 1.225 kg/m3 and kinematic viscosity 1.789410-5 kg/m2.

The equations governing the control volume are discretized. Integration of these discrete equations gives a set of equations in algebraic form. The simulations were performed using Ansys Fluent 18 R2 Computational Fluid Dynamics. The SIMPLE algorithm was chosen to correlate velocity and pressure corrections. For higher accuracy, second-order schemes of momentum, turbulent kinetic energy, and specific dissipation factor are used, and a fourth-order magnitude convergence criterion is adopted.

Table 1: Boundary conditions of problem solving

Solution methods

Fluent 18 is used in this study and the convergence criterion is 6-10. The solution algorithm is shown in Figure 6.

Figure 4: Solving algorithm

Model of Air Flow in Nasal Cavity

The most successful way to discretize the flow problem is to use the control volume method, also known as the finite volume method or Euler method. As the name suggests, the problem is split into a control volume and the conservation laws of physics apply to this volume. The governing equation is:

Continuity equation (mass conservation);

Momentum equation (Newton's second law);

Reynolds’ transport theorem.

Reynolds' transport theorem is used to convert the governing equations to Euler form [37].

Conservation of mass means that the amount of mass remains constant over time. Therefore, if there is mass flow through a volume, the accumulation rate within that volume is equal to the net outflow. In vector notation:

(1)

This is also called the continuity equation for compressible fluids in unsteady flow. For an incompressible fluid, the continuity equation is:

(2)

Reynolds’ Transport Theory

The rate of change of the general property per unit volume can be expressed similarly to the rate of change of mass per unit volume (Equation (1)) as:

(3)

The characteristic rate of change is equal to the sum of the control volume rate of increase and the net outflow.

The Momentum Equation

Once Reynolds’ transport theory is established, it can be applied to Newton's second law to compute the rate of change of momentum using a control volume approach: 

(4)

Using:

	(5)
	(6)

The total force acting on a volume is the sum of pressure, gravity (and other body forces), shear force and normal stress [37]. Applying Reynolds' transport theorem and Newton's second law to a three-dimensional infinitesimal control volume yields the three-dimensional Navier–Stokes differential equation for laminar flow [38]:

)

(7)

Due to the turbulence in this case, velocity fluctuations and additional shear stresses, the so-called Reynolds stresses, must be taken into account. Velocity is modeled as the sum of mean and variation, u = U + u’. This produces a slightly different version of Navier–Stokes:

(8)

Using the Boussinesq approximation for the Reynolds stress term (turbulent stress) on the right-hand side produces:

(9)

After insertion and rearrangement, the Reynolds-averaged Navier–Stokes equations for turbulence are obtained [39].

For steady state, the transient term on the left is of course zero, so the convective term on the left is left alone. P is the pressure; k is the kinetic energy and the remaining two terms are diffusion (turbulent stress) and viscous stress. is the Kronecker delta required for the formula to give correct results for normal Reynolds stress. 

The turbulent kinetic energy equation is:

(10)

The pseudo-vortices equation () is:

(11)

Note that the addition notation i, j = 1, 2, 3 denotes the components of the velocity vector and spatial coordinates. Turbulent viscosities and their formulations are described in detail in ref. [39]. Note that the first terms in Equations (6)– (8) are ignored since steady-state conditions are applied in this study.

Numerical Solution

Combining the governing equations yields a closed system with an equal number of equations and unknowns. What now remains is a way to solve these equations numerically. The Ansys Fluent 18 solver for incompressible steady-state turbulence, called Simple Algorithm, uses the SIMPLE algorithm to solve the flow. 

The Simple Method

SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) starts with a (preferably qualified) guess of initial and boundary conditions for the pressure field. The guessed value is marked with an asterisk, . is used to the discretized momentum equations to calculate velocities, . These velocities are referred to and treated as guessed values.

The continuity equation is used to express a correction to the pressure, . The guessed and correction values are summed to get the actual pressure. Similarly, a velocity correction is added to :

	(12)
	(13)

After correcting for pressure and velocity, the other discretized transport equations are solved. This completes one iteration. The process is repeated using the corrected pressures and velocities as initial estimates for the next iteration. This continues until the convergence criterion is met. Different criteria may exist for different flow characteristics. The convergence criteria for this paper are described in Section 2.3.7. A flowchart of the SIMPLE method by Versteeg and Malalasekera is shown below [40].

There are several versions of this algorithm, for example, SIMPLE revised or the SIMPLER method, which directly computes the pressure field via the pressure discretization equation [40]. Nevertheless, only the original SIMPLE method was used in this study, as it forms the basis for the chosen software solver Simple.

Under-Relaxation

The SIMPLE method is prone to instability. One way to prevent this is by introducing under-relaxation coefficients,, for the pressure and velocity corrections. An under-relaxation coefficient is a factor between 0 and 1 that is used to dampen corrected values:

	(14)
	(15)

Even if the solution oscillates, the amplitude will be small, so the solution may not oscillate or diverge. However, this slows down convergence. The optimal relaxation factor depends on mesh, flow problem and iteration method and is difficult to predict [40].

Numerical Schemes

The most important basic properties of discretization methods are:

Conservatism: The same number of properties entering the cell area on one side exit the same area on the other side;

Boundedness: The property value of a point is within the range spanned by the boundary

Boundary   < value At Node < boundary

(16)

Transportiveness: The influence of the property ∅ on the value of a cell by neighboring cells depends on the convective–diffusion relationship in the flow. This relationship is called the Peclet number, Pe. The influence of neighboring cells should be influenced according to the Peclet number.

(17)

For pure diffusion, −Pe > 0; for pure convection, −Pe > 1;

To use routine comparison;

The purchase number scheme is similar to other purchases.

It is hard to complete all the processes using one product. The desired properties, such as second-order precision and boundedness, can be mutually exclusive. Quality also comes at a cost. Sometimes the best systems are too expensive. However, the scheme works well, even if it has only first-order precision or is less portable, as long as it is suitable for the problem at hand. For example, the central derivative scheme is sufficient for pure diffusion, but becomes less transportable as increases [40]. In this case, the scheme was chosen to achieve second-order accuracy. 

Boundary Conditions

To create a physical or computational model, you must identify limits and set appropriate values for those limits. The nasal airways have boundaries that change over time and are usually lined with a thin layer of mucus. The shape of the nose also affects the internal flow due to the inspiratory cycle and the flow entering the nostrils. During the respiratory cycle, the airway is imprinted under the nasopharynx. Due to the aforementioned complexity and the authors' knowledge, no model yet covers all these features. In the current study, three different constant inlet velocities of 15, 17.4 and 20 LPM were used as inlet conditions to simulate normal breathing in sleep, rest and relaxation situations. In addition, atmospheric pressure conditions are set at the outlet. As already mentioned, the effects of gravitational acceleration are ignored and a no-slip condition is used on the walls of the model. 

The initial values for k and ω at the inlet are determined using the following relations [41]:

	(18)
	(19)

where is the hydraulic diameter and I is the turbulence intensity obtained from the following equation [40]:

(20)

The nasal wall is considered rigid and non-slip [31,42]. A mass flow rate is applied to the nostril inlet, which corresponds to a specific air inlet (liters/minute). The flow range considered in this study is 15, 17.4 and 20 LPM and exists in the laminar and turbulent regimes for both cases [9,43,44]. As a general rule, for an adult nasal cavity, flow below 15 LPM is considered laminar, between 15 and 17.4 LPM is transient and above 17.4 LPM is turbulent. For the nasopharynx, the boundary condition “outflow” is considered, which assumes a fully developed flow. This simulation does not take into account the presence of mucus layers and nose hairs. The physical characteristics of the fluid (air) used in this calculation are a density of 1.225 and a kinematic viscosity of 1.7894 × .

The governing equations are discretized over the control volume. Integrating these discretized equations yields a set of equations in algebraic form. Simulations were performed using the CFD solver Ansys Fluent 18 R2. The SIMPLE algorithm was chosen to correlate velocity and pressure corrections. For higher accuracy, second-order schemes of momentum, turbulent kinetic energy and specific dissipation factor are used, and a convergence criterion of four orders of magnitude is adopted.

Solution Methods

In this study, Fluent 18 is applied and the convergence criterion is . The solution algorithms are shown in Figure 6.

The initial values ​​of k and ω in the input are determined using the following relationships [41].

	(18)
	(19)

Where is the hydraulic diameter and I is the turbulence intensity, which is obtained from the following equation.

(20)

The nasal wall is considered rigid and non-slip [30 and 42]. The mass flow rate applied to the nasal cavity inlet corresponds to a specific air inlet (liters per minute or LPM). The range of flow rate considered in this study is 15, 17.4 and 20 LPM and it exists in transient and turbulent regime for both cases [9-43-44]. As a general rule, for the nasal cavity of adults, a flow of less than 15 LPM is considered calm, between 15 and 17.4 LPM transient, and above 17.4 LPM turbulent. For the exit of the nasal cavity (nasopharynx), the "exit" boundary condition is considered, assuming a fully developed flow. In this simulation, we do not consider the existence of mucous layers and nasal villi. The physical characteristics of the fluid (air) used in this calculation are density 1.225 kg/m3 and kinematic viscosity 1.789410-5 kg/m2.

The equations governing the control volume are discretized. Integration of these discrete equations gives a set of equations in algebraic form. The simulations were performed using Ansys Fluent 18 R2 Computational Fluid Dynamics. The SIMPLE algorithm was chosen to correlate velocity and pressure corrections. For higher accuracy, second-order schemes of momentum, turbulent kinetic energy, and specific dissipation factor are used, and a fourth-order magnitude convergence criterion is adopted.

Table 1: Boundary conditions of problem solving

Solution methods

Fluent 18 is used in this study and the convergence criterion is 6-10. The solution algorithm is shown in Figure 6.

Figure 4: Solving algorithm

In this research, An Iranian patient with nasal septal deviation and right nasal cavity obstruction was evaluated. The patient did not have facial trauma, respiratory infection, polyps or seasonal allergies within six months before surgery. Septoplasty was done. CT scan and 3D reconstruction of the nasal cavity model performed before operation and also three months after septoplasty.

Figure 5 shows the changes in the cross-sectional area of the nasal cavity from anterior to posterior before and after surgery. We also compared our results with previous researchers for a more detailed analysis. The cross-sectional area of ​​the nasal inlet in right nasal cavity for the pre-operative model was 1.6571 cm2. After surgery, the cross-sectional did not change at inlet.

Figure 5: Comparison of the cross-sectional area of ​​the patient's nasal cavity before and after septoplasty

Also, the cross-sectional area of ​​the nasal cavity did not change between 6 cm from inlet to nasopharynx after surgery, because this area did not have effect on obstructing the airflow according to CFD data. Increasing the cross-sectional area was at the 3 to 6 cm from nasal tip after septoplasty where this area was turbulent in CFD analysis.

Velocity​​distribution and airflow lines

Abnormal distribution and irregular patterns of the airflow lines of the nasal passages may disturb the stimulation of the olfactory area and cause 

excessive dryness of the nose and also difficulty in breathing. Figure (6) shows the paths of inhaled air for the flow rate of 17.4 LPM and it is colored according to the value of the velocity, which is displayed from the blue spectrum at the value of 1.1 m/s to the red spectrum at the value of 10 m/s. This image shows the difference in breathing flow patterns before and after septoplasty. As can be seen, in the post-operative model a larger amount of airflow passes through the middle and upper parts of nasal cavity compared to pre- operative model. Also, the narrowed area in the lower area diverts the air towards the middle and upper flow. However, the flow film within the nasal passage converges significantly after septoplasty.

Figure 6: Air flow lines at a flow rate of 17.4 liters per minute (a) before surgery (b) after surgery

In order to better understand the airflow in different parts of nasal cavity and to compare the changes in different cross-sectional areas, we considered cross-sections at approximately 1 cm intervals. Finally, we extracted the characteristics and behavior of the flow as it is shown in Figures 7 and 8.

Figure 7: Cross-sectional cuts along the nasal cavity before surgery

Figure 8: Velocity contours at a flow rate of 17.4 LPM for the nasal cavity before surgery

After the numerical analysis with the help of CFD and the diagnosis of critical points, the surgeon was asked to perform the surgery to remove the obstruction and correct the nasal septum based on the critical points of the numerical model CFD. After the patient's recovery (3 months), the steps of tomography imaging of the patient's face, 3D formation of the patient's nasal 

cavity and finally numerical analysis were performed on the modified nasal cavity of the patient and the following results were obtained.

Figure 9 shows seven coronal sections after surgery. The airflow lines are vertical to the coronal sections of the nasal cavity.

Figure 9: Coronal sections of the nasal cavity (after surgery)

Considering reconstructing the nasal cavity model after the operation, uniformity of the airflow velocity was achieved. In the distance of 3 to 6 cm from nasal tip, the velocity decreases with along the coronal sections of the right nasal cavity (figure 10).

Figure 10: Velocity contours at a flow rate of 17.4 LPM minute for the nasal cavity after surgery

Pressure profiles

When the airflow reaches the deviated septum at a constant rate in the nasal cavity, the velocity increases due to the reduction of the cross-sectional area. As a result of this behavior, the fluid pressure and static pressure decrease. Thus, the pressure difference in the obstruction area increases. Because of the flow velocity and pressure behavior in the nasal septal deviated area, the adhesion of the flow layers near the wall and creation of small rotation 

vectors in the area under the viscous layer leads to increase the shear stress and airflow regime turbulence. 

Figure 11 shows the static pressure distribution in the nasal cavity at the flow rate of 17.4 and 25 LPM. Static pressure decreases in the 3 to 6 cm from nasal tip, where the nasal septum on the right side is deviated. Value of the pressure increases at posterior sections due to the connection of the two nasal cavities.

Figure 11: Static pressure in the nasal cavity model before and after surgery

Also, the pressure drop is greater in right nasal cavity than in left cavity, which causes greater resistance to the airflow on the right side. Intranasal pressure is a measure to report nasal function. As expected, and shown in Figure 11, by correcting the obstruction of nasal cavity with surgery, the pressure drops decreases and this indicates that the person has less effort to breathe. Also, by increasing the airflow rate, the pressure drops increases along the coronal sections before and after the surgery.

Figure 11: Pressure difference in the nasal cavity model before and after surgery

Wall Shear Stress (WSS)

The shear stress of the wall along the nasal cavity means the adhesion of the airflow layers near the surface to the wall, When the velocity increases according to the equation of shear stress (21), the shear stress also increases. As figure 12 shows WSS increases at the level of inferior turbinate in both before and after surgery models.

(21)

WSS for the nasal cavity model with NSD had the highest value of 0.04 Pa at 4 cm from nasal tip. After septoplasty, the highest WSS was at 5 cm from nasal tip and decreased 0.0278 Pa in the mass flow rate of 17.4 LPM (Figure 12 shows WSS in both pre and post septoplasty models)

Figure 12: Average wall shear stress distribution before and after surgery

Airflow resistance

To calculate the value of resistance (R), we used Zhu et al.'s formula [45 and 46].

  (22)

where ∆P is the pressure difference between the nasal cavity and the desired cross section and Q is the flow rate in the nasal cavity. According to the above equation, the resistance in the models before and after surgery at the flow rate of 17.4 LPM was 0.977 and 0.49 respectively. Figure 13 shows resistance in different mass flow rates.

Figure 13: Air flow resistance in the right nasal cavity at different flow rates

Hydraulic power

The hydraulic power in the nasal cavity was evaluated because decreasing of this parameter leads to regime of flow changes from turbulence to laminar. In the nasal cavity which had a nasal deviation, the hydraulic power increased in the area of obstruction. 

Hydraulic power can be defined as follows:

W= ∆pQ (23)

Figure 14 shows the hydraulic power before and after surgery at the level of obstruction. As can be seen, due to the increase in the cross-sectional area of ​​the right nasal cavity and the removal of the obstruction after septoplasty, the hydraulic power at different flow rates decreases, which indicates decreasing the patient's breathing effort.

Figure 14: Hydraulic power of nasal cavity model before and after surgery

It was observed that by using CFD before surgery and consultation with the surgeon based on the obtained results, the healthy nasal tissues of the patient, which do not need to be destroyed during surgery, remain anatomically, which indicates the usefulness of this objective method.

Deviation of septum is one of the most common causes of nasal airway obstruction. Nasal septal deviation plays an important role in nasal anatomy [46]. The anatomy of the nasal cavity, which is a conical channel with three lower, middle and upper turbinates on the lateral wall is complex. For this reason, many researchers use objective methods to examine the upper airway and plan surgery. Hariri et al studied the relationship between CFD and rhinomanometric results in septoplasty. They found that the results of pressure drop obtained from CFD are well correlated with the results of preoperative rhinomanometry, 

CFD may be able to detect specific obstructed sites and post-surgical results [47]. Dmitry Teri Tyakov and colleagues discussed potential applications of CFD in otolaryngology, including assessment of nasal airway abnormalities and surgical interventions. They also pointed to limitations in current modeling approaches [34]. 

A study by Springett et al on how CFD modeling with 3D printing and virtual anatomy tools has the potential to make surgical planning is done. They showed that by using CFD and 3D printing highly personalized surgery can be planned [48]. 

One of the diseases that affects the nasal cavity is the perforation of the nasal septum, which has various causes, mainly trauma or systemic diseases. Considering that the methods and effects of different treatments are challenges for specialists, some researchers such as Won et al. and Kimble et al. used CFD models to answer questions related to this pathology [49, 51]. 

Treatment of nasal obstruction due to deviated nasal septum is another challenge that CFD can help to choose the best treatment. For example, nasal turbinoplasty or lower turbinate size reduction is a procedure that is sometimes performed on patients with nasal obstruction. The purpose of reducing the turbinates is to improve nasal airflow by reducing the size of the concha. CFD studies have shown that nasal turbinoplasty can increase nasal airflow, reduce nasal resistance and improve temperature and humidity regulation in the nasal cavity [35, 51].

It is useful to mention that the modeling based on CFD often considers the flow to be incompressible and the walls to be rigid, which are actually inherent uncertainties. Also, because of differences in methods and sample characteristics, comparisons with normative values ​​should be interpreted with caution. According to the above points, in our study, the velocity distribution pattern in the patient with deviated nasal septum shows an unusual velocity, which is different from the normal nasal cavity, where the maximum velocity is located in the area of ​​the nasal valve [47]. In this study, the location of the maximum velocity was changed to the nasal valve after septoplasty. This fact shows that septoplasty based on CFD data can normalize the velocity distribution. After surgery, the area of ​​the nasal valve becomes the narrowest part of the nasal cavity, as a result of which the speed along the nasal cavity decreases. According to Croce et al.'s study and our study, there is a sudden increase in cross-sectional area after the nasal valve, which leads to a decrease in velocity. Nasal cavity with deviated septum experiences maximum velocity in the mid-nasal region, which indicates an abnormal flow pattern. This can potentially affect the thermal ventilation and filtration capabilities of the nasal cavity [52].

Pressure drop is traditionally used to validate CFD solutions, measuring the pressure difference between the nostril and the nasopharyngeal exit [50]. According to a study by Ottaviano and Fokkens, pressure drop is thought to have a significant impact on nasal patency [53]. 

In this study, we observed higher values of pressure drop in pre operative nasal cavity with septal deviation compared to post operative. This issue is in agreement with previous studies [49, 53].

J. Corda, et al. reported higher pressure drop and airflow disturbance in deviated nasal cavity and also noted the effectiveness of using CFD modeling to show obstruction caused by deviated nasal septum in these patients [54].

WSS is the resistance force that is produced when the air flows within the nasal passage, allowing for the exchange of heat and moisture between the air and the nasal pasage. [22].

In this study, we showed that by specifying the critical points of the deviated nasal cavity and suggesting to the surgeon to do septoplasty according to CFD data, the WSS reduced. 

WSS and the location where the fluid has the highest adhesion have been well discussed in Corda et al.'s study [53]. Their results showed that the nasal valve has the highest wall shear stress. Overall, our findings are in good agreement with their results, in a way that most researchers share the same opinion.

A study by Croce A. et al. compares WSS, velocity profiles and pressure drop in NSD patient pre and post septoplasty. All of mentioned data had significant difference after septoplasty. This finding aligns with the present study’s results [54].

If we consider the function of airflow passing through the nasal cavity after septoplasty surgery to be similar to that of a healthy nasal cavity, in the model after septoplasty, a boundary layer is formed like a healthy nasal cavity, but in the obstructed nasal cavity, the boundary layer does not produce due to the sudden separation of the flow. For this reason, we used k-ω SST Transient (y+) method in solving the problem of this research.

Our study showed relatively higher resistance values ​​in the preoperative nasal cavity with NSD compared to the postoperative nasal cavity.

The resistance in the nasal cavity of a healthy person has been studied in previous studies. Wen et al. [22], Zubair et al. [55], and Weinhold and Mlynski [56] reported the normal resistance of the nasal cavity at a flow rate of 20 LPM as 0.054, 0.06, and 0.068, respectively [22, 55, 56].

The sense of obstruction can be clarified in the case of NSD because the resistance is greater than in the normal nasal cavity.

While many studies report higher resistance in cases of airway obstruction, not all do. Further research is needed to clarify which patient and septum specific factors most influence on resistance.

Deviation of the septum as one of the main causes of nasal obstruction, was evaluated in a patient with NSD before and after septoplasty. In this study, pre and post operative CT scan of was used to make a 3D model and airflow analysis at 15, 17.4 and 20 LPM was performed. The narrowest area in the case of nasal septum deviation was at the middle of the nasal cavity before septoplasty, while this area changed to the entrance (valve) of the nose after septoplasty. Velocity patterns showed the maximum in the middle nasal cavity in case of preoperative nasal model, which may affect the physiological functions of the nose. The model before surgery shows higher resistance than the nasal cavity after septoplasty. The values ​​of pressure drop in the nasal cavity before the operation were higher than the nasal cavity after the operation. Also, WSS and the hydraulic power were reduced after septoplasty. 

Therefore, using CFD data for determining the exact area that caused obstruction and suggesting the surgeon to do septoplasty according to mentioned area not only saves the healthy nasal tissue, but also improves the results of septoplasty. 

Institutional Review Board Statement: The study was conducted in accordance with the Declara-tion of Mashhad university of medical sciences, Iran, and approved by the Institutional Review Board

Informed Consent Statement: Informed consent was obtained from subject involved in the study.

Data Availability Statement: Data is contained within the article.

Deviation of the septum as one of the main causes of nasal obstruction, was evaluated in a patient with NSD before and after septoplasty. In this study, pre and post operative CT scan of was used to make a 3D model and airflow analysis at 15, 17.4 and 20 LPM was performed. The narrowest area in the case of nasal septum deviation was at the middle of the nasal cavity before septoplasty, while this area changed to the entrance (valve) of the nose after septoplasty. Velocity patterns showed the maximum in the middle nasal cavity in case of preoperative nasal model, which may affect the physiological functions of the nose. The model before surgery shows higher resistance than the nasal cavity after septoplasty. The values ​​of pressure drop in the nasal cavity before the operation were higher than the nasal cavity after the operation. Also, WSS and the hydraulic power were reduced after septoplasty. 

Therefore, using CFD data for determining the exact area that caused obstruction and suggesting the surgeon to do septoplasty according to mentioned area not only saves the healthy nasal tissue, but also improves the results of septoplasty. 

Institutional Review Board Statement: The study was conducted in accordance with the Declara-tion of Mashhad university of medical sciences, Iran, and approved by the Institutional Review Board

Informed Consent Statement: Informed consent was obtained from subject involved in the study.

Data Availability Statement: Data is contained within the article.

The authors declare no conflict of interest

The authors declare no conflict of interest