OPEN-SOURCE CFD ANALYSIS OF NASAL FLOWS

. Recently, as studying the transmission of airborne particles and the effect of filtration devices became of pronounced importance due to the COVID-19 pandemic, having the ability to visualise nasal flows and perform numerical analyses may prove to be an invaluable resource. What is more, due to the complexity of the anatomy of the nasal cavity and its variations in patients, a patient-specific analysis could help in treating different conditions affecting nasal flows. That is why this research describes a democratised and easily reproducible procedure used as a tool for the investigation of nasal flows using open-source computational fluid dynamics. A detailed overview of a procedure developed for analysis of nasal flows is presented. Every step of the procedure, from creating the computational model to airflow simulation, is based on freely available and open-source tools, thus providing a widely available framework for investigation of nasal flows. The framework consists of extracting a computational model from computed tomography (CT) or magnetic resonance imaging (MRI) scans, model refinement and numerical simulations performed using foam-extend . The framework was applied to real radiological (CT) scans of the nasal cavity and was used for numerical simulations of a transient, incompressible, turbulent, single-phase nasal airflow. The simulations were performed using a novel approach based on the pUCoupledFoam coupled incompressible pressure-velocity solver available in foam-extend , used in conjunction with Block-SAMG - a block-selective algebraic multigrid solution algorithm.The analysis of the results of transient simulations provided insight into the behaviour of the nasal airflow in a real (patient-specific) geometry. Generally, the simulations show how the flow is affected by anatomical barriers, resulting in higher velocities throughout the left side of the particular anatomy in question. Furthermore, the analysed cross-sections showed the greatest values of the velocity magnitudes appearing along the nasal septum (nasal wall separating the cavity), with significantly smaller values noticeable in the peripheral region of the geometry.


Introduction
The recent global pandemic prompted the scientific community to focus their attention to the investigation of the transmission of airborne particles, testing of existing and the development of improved personal protection equipment (PPE). A key factor in both the investigation of particle transmission and PPE is human breathing, more precisely nasal and oral airflows.
Breathing is one of the primary functions of the human body. The nose and the nasal cavity are part of the upper respiratory tract and, as such, are indispensable for the breathing process. The nasal cavity has a complex role in the breathing process, taking part in temperature and humidity regulation of the inhaled air and the filtration of particles. Not only is the geometry of the nasal cavity and the surrounding tissue highly complex, but the geometry differs widely for each person. Being able to accurately determine the specifics of airflow in the complex geometry of the nasal cavity, would provide numerous insights into phenomena associated with the breathing process, such as the effects of nasal airway obstructions, transport of particles and the investigation of recirculation zones. Each of the aforementioned phenomena can have an adverse effect on the breathing process, the health and comfort of a patient and the choice of treatment.
Developments in the field of radiology and medical diagnostics led to the possibility of creating threedimensional models of the nasal cavity without the need for invasive procedures. Furthermore, advancements in the field of computational fluid dynamics (CFD) made it possible to achieve highly accurate numerical analyses of laminar and turbulent flows in complex geometries [1]. This makes it possible to use CFD as a tool for investigating the complexities of nasal airflows.
Early attempts at using numerical analysis for nasal flows started in the 90s. Due to limited computer resources, the computational models used were simplified [2] or only part of the nasal cavity was considered [3]. With the development of computational resources, more researchers begin using numerical analysis and CFD to investigate nasal flows [4,5,6,7,8,9]. More recently, CFD was used to investigate different aspects of nasal flow: particle deposition in steady and pulsating flows [10,11], improving nasal surgery results by investigating the effects of surgical procedures [12], Large Eddy Simulation (LES) of water droplet deposition and assessment of nasal obstructions [13]. While most of the aforementioned studies relied on proprietary software for meshing and numerical analysis, Tretiakow et al. [14] used OpenFOAM to perform a series of CFD simulations of airflow in the nasal cavity and paranasal sinuses. Tretiakow et al. used a workflow based on open-source and freeware software, but also used an external meshing software, which requires additional steps in the simulation procedure.
In this paper, a detailed overview of a procedure developed for analysis of nasal flows is presented. The entire procedure, from the creation of a computational model to the airflow simulation, is based on opensource software, providing a widely available tool for detailed investigation of nasal flows. For the first time a combination of the pUCoupledFoam coupled incompressible pressure-velocity solver in conjunction with the Block-SAMG block-selective algebraic multigrid solution algorithm is applied to biological flows, resulting in a novel approach to simulating nasal flows. This study was performed in hopes of creating a democratised procedure based on non-commercial software which would allow for a wider application of CFD simulations of nasal airflow in the investigation of PPE, particle transmission and the in application of patient-specific medical procedures, including virtual surgery.

Workflow Overview
A simplified workflow used for conducting a single nasal airflow simulation is outlined by the steps given here: • Extracting a 3D model of the nasal cavity from a CT/MRI scan. • Refining the model as to obtain a airtight geometry suitable for the meshing algorithm.
• Creating an unstructured computational mesh using an automatic meshing algorithm.
• Setting up the appropriate boundary conditions.
• Running the simulation in foam-extend.
• Processing the results: visualisation and analysis. A more detailed overview, providing further information and comments for each step outlined here, is given in Appendix B.
The initial step of the procedure requires a computed tomography (CT) or magnetic resonance imaging (MRI) scan of the anatomy of interest, i.e. the nasal cavity and the surrounding tissue. The radiological data from a CT/MRI scan may then be processed in 3D Slicer -a medical image analysis and volume rendering platform [15]. Once exported in a suitable format, the rendered volume is used to create a computational domain with the help of a meshing software. In this workflow, the computational domain was generated using cfMesh -a cross-platform library for automatic mesh generation [16]. The final steps of the workflow focus on creating a computational domain, setting up numerical simulations for different boundary conditions and analysing the results.

Model Extraction and Refinement
Radiological data (e.g. a CT scan) depicts the scanned tissue as a series of axial slices composed of thousands of pixels (square elements with an area of around 1 mm 2 ), where each pixel is assigned a greyscale value between 0 and 255 (from black to white) [17]. Furthermore, the axial spacing between each slice usually equals 1 mm, thus transforming each pixel into a three-dimensional element, a cube-shaped volume -a voxel. Voxels that represent different tissue have different values of the Hounsfield Unit (HU). The authors in [18] provide a detailed overview of how different voxels are assigned different HU values depending on the levels of x-ray attenuation. Knowing the HU values for different tissue allows us to separate and extract the voxels from a set of two-dimensional slices, which represent a three-dimensional area of interest. This process is called volume-rendering. An example of the volume rendering process in 3D Slicer, depicting three anatomical planes and the rendered volume, is shown in Fig. 1. The rendered volume shown in Fig. 1 (green region) represents the volume filled by air in the region of interest, i.e. the air inside the nasal cavity and in front of the outer nose.
The volume-rendering part of the workflow is carried out in 3D Slicer [15]. This open-source utility was used to transform a series of two-dimensional CT scans (slices) into a three-dimensional representation of the nasal cavity. The volume-rendering procedure requires CT scans of the area of interest to be taken in the transverse, frontal and sagittal anatomical planes. 3D Slicer allows for the CT scans situated in the aforementioned planes to be superposed and the resulting three-dimensional representation of the tissue to be rendered as a volume. As different types of tissue correspond to different values on the HU scale, different parts of the anatomy may be selected for rendering by changing the HU range used during rendering. The final rendered volume can be exported as an STL file. Before the rendered STL model can be used in in the mesh generation procedure, the surface must be edited using a tool capable of editing STL surfaces. For this project, the surface was modified in Blender -an open-source 3D computer graphics toolset [19]. Blender was used to remove unwanted parts of the anatomy that could not be excluded during the volume rendering procedure and to make the STL model airtight. Furthermore, a spherical inlet patch was added to represent the atmosphere and the STL surface was divided into separate patches (need for defining the boundary conditions). The final STL used in the meshing procedure is shown in Fig. 2.

Computational Domain
The spatial discretisation (meshing) was done in an unstructured manner using the cartesianMesh algorithm in cfMesh [16]. Further information regarding the meshing procedure may be found in the Appendix B.
The size of computational domain was determined so that the mesh resolution used is fine enough to accurately describe the features of the nasal airflow and with sufficient cell resolution even in the more complex regions of the geometry, which are characterised by very small cross-sections. On the other hand, mesh resolution was limited to allow the workflow to be used without the need for high-performance computing (HPC) equipment, while maintaining acceptable time-to-solution.
The resulting computational domain is an unstructured, predominantly hexahedral mesh, consisting of 4 407 571 control volumes. The maximum control volume size is roughly 1 mm, while the boundary cell size near the walls of the nasal cavity is approximately 0.3 mm. A depiction of the final computational  mesh, halved by the sagittal (longitudinal) plane, is depicted in Fig. 3, while a detailed overview of different cell types forming the mesh is given in Table 1. It should be noted that, while the computational mesh depicted in Fig. 3 may seem very structured, it is indeed an unstructured hexahedral mesh which uses polyhedra for transitioning between different cell sizes. The predominant use of hexahedral cells gives the mesh its ordered appearance. It should be noted that additional refinement was applied to the surface walls of the nasal cavity and the spherical patch representing the atmosphere to better capture the curvature of the surfaces and the transition between the outer surface of the nose and the spherical patch. While no additional boundary layers were created, care was taken during the meshing procedure to ensure that the values of y+ are within the acceptable range for use with wall functions, i.e. between 30 to 300. The readers should also note that the geometry and the resulting computational domain presented in this section could not be made freely available to the readers. That is why a simplified geometry of the nasal cavity was created and was made available. Further information regarding the simplified geometry may be found in Appendix C.  with ω as reference value

Numerical Simulation
The previously described computational domain was used as basis for running three distinctive sets of simulations: an inhalation lasting 2 seconds, an exhalation lasting 3 seconds and a quick inhalation ("sniff") lasting 0.5 seconds. All simulations were performed using foam-extend [1], a community driven fork of the OpenFOAM library for general computational continuum mechanics [20].
The nasal airflow was modelled as a transient, incompressible, turbulent, single-phase flow. The simulations were run using the implicitly-coupled incompressible pressure-velocity algorithm described in [21]. Previous research showed that the implicitly-coupled procedure is more stable for cases with internal fluid flow and that the time-to-solution is significantly reduced in comparison to segregated algorithms, such as the SIMPLE method. Details of the Block-SAMG algorithm are given in [21,22].
Here, we shall briefly describe the recommender settings for optimal convergence of the block pressurevelocity linearised system for future reference and as best practice advice for foam-extend users. Since the pressure-velocity matrix is asymmetric and consists of equations of a different type, namely hyperbolic for the momentum equation and elliptic for the pressure equation, the choice for the linear algorithm is not straightforward. With segregated methods, the pressure equation is solved using linear algorithms for symmetric matrices, usually conjugate gradients or multigrid methods with fixed-point smoothers (Jacobi, Gauss-Seidel), while the asymmetric momentum equation is solved with preconditioned Krylov subspace methods (GMRES or BiCGStab). The novel Block-SAMG enables an efficient solution of the pressure-velocity system, where the coarse-level creation in multigrid is algebraic, governed by the pressure equation. ILUC0, a typical preconditioning method for the momentum equation, is used as a smoother with Block-SAMG, even though its smoothing properties are not mathematically proven in literature. Thus, the proof is numerical and is presented in [22,21]. Turbulence was modelled using the standard k − ω − SST model available in foam-extend as the model proved to be stable, allows for modelling turbulence near the boundary and the flow far from the wall. The model produces very good results for incompressible internal and external flows. The k − ω − SST is implemented in foam-extend according to the model presented by Menter and Esch [23]. Furthermore, the update coefficients for the model are used according to [14], while the consistent production terms from [23,24] are updated according to [25].
All three simulations were set up with a time-step of 0.001 s, which may be seen in the provided case files. The resulting execution times for each simulation using the aforementioned value of the time-step and running are presented in Table 3. The value of 0.001 s as the time-step results in the mean value of the Courant number of approximately 0.58 for inhalation, 0.5 for exhalation and 1.05 for the short inhalation ("sniff"), while the value of the local Courant number varies at each time-step. The values of the Courant number were provided for informative purposes as we are not limited by the Courant number due to the fact that the simulations are run using a block-coupled solver.

Boundary Conditions.
Due to the inability to correctly asses flow conditions at each nostril, a spherical boundary was added (shown as Atmosphere in Fig. 2), with conditions corresponding to normal atmospheric conditions. Time-varying values of velocity were used at the outlet of the nasal cavity (shown as Trachea in Fig. 2). The complete list of boundary conditions (BCs) is given in Table 2.
A time-varying boundary condition for airflow velocity was used on the outlet patch (trachea) of the nasal cavity. The time-varying values of the outlet velocity were interpolated from a dataset based on different breathing curves, representing airflow through a normal (healthy) nasal cavity. Depending on the transient simulation (inhalation, exhalation or quick inhalation) a different dataset was used.
(timeVaryingUniformFixedValue) -a standard boundary condition available in foam-extend was used to facilitate reading and interpolation of the outlet velocity data. The aforementioned boundary condition is a time-varying form of a uniform fixed value BC which uses interpolation tables to interpolate the value to be applied from a supplied data series. The applied outlet boundary condition specifies a fixed value across the outlet patch at each time-step, thus imposing plug flow conditions at the outlet.
The time-varying values of velocity used at the outlet patch were derived from breathing curves depicting volumetric flow rates through the nasal cavity. Breathing curves are usually determined experimentally for a specific patient using different clinical measurement techniques such as acoustic rhinometry or hot-wire anemometry [27]. Due to the fact that the authors of this paper were unable to obtain the exact breathing curve (nasal flow measurements ) of the patient whose radiological data was used to create the computational model of the nasal cavity, data available in the literature was used. Furthermore, breathing curves presented in various literature usually depict values of the volumetric flow rate for a specified duration of inhalation or exhalation and derived from measurements at the nostrils (inlet). As incompressible flow is assumed, the same volumetric flow rate could be assumed at the cross-section at which the outlet BC was placed (trachea) and, knowing the area of the outlet cross-section, the velocity at the outlet boundary could be deduced.
Volumetric flow rates (breathing curves) used to calculate airflow velocities through the outlet patch were found in literature. More precisely, the data for a quick inhalation ("sniff") was presented by Bates et al. [30], who used a breathing curve derived from the experimental data of Rennie et al. [27]. Experimental measurements carried out by Rennie et al. used different clinical techniques, such as bilateral hot-wire anemometry and acoustic rhinometry, to determine time-varying profiles of nasal inspiration at rest and for "sniffing". The authors from [29] used breathing curves to set the volumetric flow through a physical model of the nasal cavity before and after a nasal turbinectomy. The breathing curves presented in [29] are referenced as experimental data found in The Handbook of Physiology [28]. As the source material could not be found, a breathing curve representing 5 second inhalation/exhalation cycle presented in [29], was used as input for the BCs in the case or normal inhalation/exhalation presented in this paper. The exact breathing curves used for the time-varying BC at the outlet are shown in Figures 4 and 5.
Each case (inhalation, exhalation and "sniff") was viewed as a separate and independent event, thus zero value of velocity was assumed throughout the domain at T = 0 s. This would correspond to conditions where breath is held by the patient before each inhalation, exhalation or "sniff". Furthermore, as data from the literature used for specifying the outlet BCs did not provide information for several complete inhalation-exhalation cycles, the inertial effects of prior several breathing cycles were not analysed as part of this study.

Results
The aforementioned computational domain and boundary conditions were used to set up three distinct sets of transient simulations: • a simulation of an inhalation based on the breathing curve of a healthy adult person and the previously described geometry, with a duration of 2 seconds , i.e. a normal inhalation, • a simulation of an exhalation using the same assumptions as above and lasting for 3 seconds, i.e. a normal exhalation, • and an intense inhalation based on the breathing curve shown in Fig. 5, with a duration of 0.5 seconds, i.e. a quick inhalation -"sniff " It should be noted that the authors used different scales for depicting values of the velocity magnitude and turbulent kinetic energy throughout the figures presented next, as the use of a consistent scales would result in diminished quality of the presented results. The reader is reminded to take note of the different scales given for the presented fields. What is more, the readers are advised to use the overview of the anatomy of the nasal cavity provided in the Appendix A as a reference for a better understanding of the model and the results presented here. Values of turbulent kinetic energy (TKE) were investigated as part of this study. Figure 8 shows values of TKE at T = 1 s for four cross-sections of the nasal cavity. Relatively higher values of TKE can be seen at the throughout the left-hand side of the cavity when compared to the right-hand side, which is attributed to the higher values of airflow and velocity through the left region of the nasal cavity, which shows near-zero values of TKE. It is interesting to note that non-zero value regions of TKE may be found at the entrance into the maxillary sinuses. This may be indicative of the existence of a low-velocity flow into the maxillary sinuses, which was not captured in Fig. 7. 6.2. Normal Exhalation. Transient simulations of a normal exhalation lasting for 3 seconds were conducted and the result are presented next. Streamlines depicting airflow during exhalation for several time steps of the simulation are shown in Fig. 9. During exhalation, the direction of the flow reverses, the air enters the nasal cavity at the pharynx and exists through the nostrils. As before, the velocity of the airflow increases as the airways become narrower near the inner nasal valve and decreases once it moves past the valve. The values of velocity increase during exhalation due to an increase in the volumetric flow rate, peaking at T = 1 s. At this time, the maximum airflow velocity of 4.51 m/s is reached near the inner nasal valve. After the simulation reaches T = 1 s, the volumetric flow rate and the airflow velocity gradually decrease throughout the nasal cavity. Similar to what was observed during inhalation, there is very little flow present in the peripheral regions of the nasal cavity. Attributed to the reverse airflow   direction during exhalation, a recirculation of airflow from the nasal cavity into the sphenoidal sinuses can be observed. Fig. 10 shows several cross-sections of the flow field inside the nasal cavity during exhalation. Higher values of velocity and airflow are localised to the middle meatus region near the nasal septum, as was the case during inhalation. On the other hand, as the flow direction is reversed, the flow has more time to develop before reaching the nasal valve region, resulting in a more uniform distribution across the cross-sections of the inferior and superior meatus region. Fig. 11 depicts turbulence, represented using TKE, for the case of normal exhalation at T = 1 s. An increase in TKE values may be noted throughout different cross-section. This increase is localised in the middle meatus region of the left-hand side of the nasal cavity, which corresponds to the region of highest airflow velocity. On the other hand the values of TKE throughout the most part of the right-hand section of the nasal cavity approach zero values. 6.3. Quick Inhalation -"Sniff ". The last set of simulations was conducted for a quick inhalation -a "sniff", based on the breathing curve presented in Fig.5 The airflow is characterised by a high volumetric flow rate, achieved in a short time period (0.5 seconds). The airflow through the nasal cavity is shown in Fig. 12. As it is defined by the same geometry, the flow looks similar in shape to the flow observed for normal inhalation, but achieves higher values of velocity due to an increased volumetric flow rate. As the flow rate increases, the maximum value of velocity (10.28 m/s) is reached at T = 0.2 s. During the same time step, average flow velocity reaches 5 m/s, which is notably higher than the maximum average velocity of 3.6 m/s, reached at T = 1 s for normal inhalation. As before, the occurrence of velocity maxima is localised to region of the internal nasal valve.
As expected, the cross-sections depicted in Fig. 13, show higher flow velocities in the middle meatus. Furthermore, a region where the airflow is absent can be observed in the middle meatus region of the right side of the nasal cavity, which is again attributed to anatomical barriers.
TKE fields for T = 0.25 s are shown in Fig. 14. An increase in values of TKE can be observed in the middle meatus region of the left side of the cavity. The presence of non-zero values of TKE my be noted at the entrance of, and inside, the maxillary sinuses. A similar phenomena was also observed for normal inhalation. On the other hand, for a quick inhalation, an increase in TKE may also be noted inside the ethmoidal sinuses.

Conclusion
This paper focused on presenting a fully open-source framework developed for the investigation of airflow phenomena in real-life anatomical models of the nasal cavity. Steps were outlined, describing how freely available volume rendering and modelling tools can be used in conjunction with foam-extend to create a numerical simulation based on a computational model of the nasal cavity and to conduct different transient simulations of nasal flows. The aforementioned framework is based on a workflow which includes: the extraction of a computational model form radiological images (CT/MRI scans) using 3D Slicer, model refinement and modification in Blender, creation of a computational domain using cfMesh and running numerical simulations using foam-extend. Simulations of a transient, incompressible, turbulent, single-phase airflow through the nasal cavity were performed using the pUCoupledFoam coupled incompressible pressure-velocity solver in foam-extend. Additionally, block-selective algebraic multigrid (Block-SAMG) was used as the solution algorithm. This is the first time the pUCoupledFoam solver and the Block-SAMG solution algorithm were applied together for use on biological flows.
A computational mesh of satisfying quality was created and the appropriate boundary conditions were applied. Velocities at the outlet were set up according to the values of the volumetric flow rate through the outlet patch, and based on breathing curves found in literature. Finally, a set of three simulations was presented using breathing curves for normal inhalation, normal exhalation and a quick inhalation ("sniff"). The results were analysed and the following conclusions were reached: • Anatomically correct three-dimensional models suitable for use in computational fluid dynamics simulations (CFD) can be created using medical data such as CT/MRI scans. • For the first time, a novel approach using the pUCoupledFoam coupled incompressible pressurevelocity solver in conjunction with the Block-SAMG was applied for the analysis of biological flows. block-selective algebraic multigrid solution algorithm was successfully applied • The results of transient simulations provided valuable insight into the behaviour of the nasal airflow in a real geometry. Generally, the simulations showed higher values of velocity to be localised to region of the inner nasal valve. Due to anatomical barriers of this particular nasal cavity, the flow was reduced in the right side of the cavity, as can be seen for different cross-sections of the cavity. The cross-sections showed increased flow and velocity in the middle meatus region of the left side of the cavity, with respect to values seen on the right side.      The presented study showed that valuable insights into phenomena connected to nasal airflows may be determined through the use of CFD simulations. Furthermore, an open-source workflow can be applied to perform CFD analysis on real patient medical data, providing patient-specific information that could be used to tailor the course of treatment or provide additional data for surgical procedures.
A future, more comprehensive study should be performed to address some concerns left unanalysed in this study. Firstly, the issue of the size of the spherical inlet patch and its proximity to the nose should be addressed. Similarly, the proximity of the outlet patch to the region of interest should also be examined. Secondly, the data used for velocity boundary conditions at the outlet should be re-evaluated. This study used general data found in the literature derived from experimental measurements performed on different patients. A future study should use clinical measurements performed on the same patient whose radiological data is used to create the model of the nasal cavity. Experimental data for several breathing cycles should be collected, as to allow for the simulation of multiple breathing cycles and the analysis of inertial effects left from previous breathing cycle. The simulations preformed in this study modelled the structure of the nasal cavity as rigid, thus a future study could investigate the effects of the deformation of the soft tissue using fluid-structure interaction simulations. Finally, effects of using a laminar-turbulent transition model should be investigated.
The workflow could also be expanded to include virtual surgery. Such a procedure would consist of editing the model of the nasal cavity to mimic any changes to the real anatomy that would happen during a real-life surgical procedure, running a series of CFD simulations on the model (before and after the changes are made) and analysing the results in hope of determining the effects of a real surgical procedure. What is more, the framework could also be used to analyse how differently sized particles are carried with the nasal flow during inhalation/exhalation and to study the effects of varying protective and filtration equipment on nasal flows.

acknowledgements
The radiological images used in this study were provided by the Sisters of Charity teaching hospital (KBC "Sestre milosrdnice") in Zagreb. Their help in finding real patient radiological data is acknowledged and greatly appreciated.
Finally, we would like to express our gratitude to TeachMeAnatomy (https://teachmeanatomy.info) for allowing us to use their anatomical illustrations in our publication. Appendix A. Anatomy of the Nasal Cavity An overview of the basic anatomy is presented in this section of the appendix. This section was included as to provide a reference for the anatomical features mentioned in the main text of the paper. Fig. 15 shows a depiction of the external nose. The external nose has a pyramidal shape and allows for air to enter and/or exit the nasal cavity. The main features of the external nose, shown in Fig. 15, are the nasal root (radix nasi ) and bridge, the nasal ridge (dorsum nasi ), the tip (apex nasi ) [31]. The elliptical openings near the tip of the nose are called the nostrils of the nasal passages (nares), while the outer walls which encircle the nostrils from each side are called nasal wings (ala nasi ) [31,32].
The nasal cavity (cavitas nasi ) is shown in Fig. 16. The nasal cavity is considered as entrance into the upper respiratory tract. The inside of the nasal cavity forms the passages which allow for air inflow or outflow during breathing. The main functions of the nasal cavity include the transport of air, filtration and temperature and humidity regulation. The nostrils are the entrance into the nasal cavity, while the main volume of the cavity is split into two roughly equal halves by a wall called the septum (septum nasi ) [31]. The expanded portion of the nasal cavity which directly follows the nostrils is called the vestibulum (vestibulum nasi ) and this section of the cavity is used for filtration of particles. The nasal valves (limen nasi ) mark the transition from the vestibule into the respiratory part of the nasal cavity. Each of the two halves of the respiratory section of the nasal cavity have three curved protrusions called turbinates (conchae nasi ) which divide the nasal cavity into the upper (meatus superior ), the middle (meatus medius) and the lower (meatus inferior ) nasal passage. The nasal passages and turbinates are shown in Fig. 17 [32]. The apertures in the posterior of the nasal cavity are called internal nostrils (choanae) and they mark the transition form the nasal cavity into the upper part of the throat (nasopharynx ) [32,33].
Generally, the name paranasal sinuses is given to air-filled volumes which are placed around the nasal cavity and connect to it by passages [32]. Their exact function is still under debate, but they are most  likely used in temperature, humidity and pressure regulation. They are also likely used to reduce the weight of the head and in voice formation [33]. The paranasal sinuses may be classified by shape, size or function. Figure 18 shows the location of the main paranasal sinuses, which include: the maxillary sinuses (sinus maxillares) -located below the cheeks, the frontal sinuses (sinus frontales) -located behind the frontal bone of the forehead, the sphenoidal sinuses (sinus sphenoidales) -located in the sphenoidal bone and the ethmoid sinuses (sinus ethmoidales) -situated in upper part of the nose (ethmoid bone) between the eyes [31,33].

Appendix B. Model Creation Workflow
In this section of the Appendix, a description of the procedure used to create the computational domain is given, starting from the initial CT scans of the nasal cavity up to the point where the computational domain may is ready to be used for simulating nasal airflow in foam-extend.
The complete procedure used to prepare the computational domain consists of four steps: volume rendering, segmentation, model creation/editing and meshing. As was already mentioned, the whole procedure is carried out using only open-source and freely available software, making it easily reproducible for anyone interested in creating a computational model from CT/MRI scans. B.1. Volume Rendering and Segmentation. The initial part of the procedure used to create a computational model from radiological data (CT/MRI scans) consists of the volume rendering procedure and segmentation of the rendered volume. Volume rendering is the process of visualising volumes as 3D object and segmentation is the process if delineating structures or interest. Both procedures are carried out using 3D Slicer -a medical image analysis and volume rendering platform [15].
The graphical user interface of 3D Slicer offers a combined view of the currently open module and four different representations of the currently loaded data, as shown in figure Fig. 19. The four available views are the axial (Fig. 19 -red window), the sagittal (Fig. 19 -yellow window), the coronal (Fig. 19 green window) and the three-dimensional view (Fig. 19 -blue window).  Fig. 19 can be used to import different datasets, including radiological data. Medical scans, such as CRT or MRI radiological data, are usually provided as DICOM data files. DICOM objects are the standard way of storing and transporting medical data, radiological scans and patient information. The Welcome module may be used to load DICOM objects containing radiological scans which are to be used in the volume rendering procedure. Fig. 20 shows an redacted example of the information contained in a DICOM data file and the CT scans provided as part of the DICOM object. Once the DICOM data file has been successfully imported, a series of radiological scans may be loaded to provide the basis for the rendering of a 3D volume.
After a series of CT/MRI scans is selected and loaded successfully, the data may be viewed in the three plains shown in Fig. 19 (axial, sagittal, coronal), but needs to be rendered first to be viewed as a 3D volume. The volume rendering procedure is procedure is performed using the Volume Rendering Module, allowing for several properties to be changed to provide the best representation of the anatomy of interest. The most important properties in the Volume Rendering Module are shown in Fig. 19 and include: Volume -allowing for the selection of the data which will be used for rendering, ROI -allowing for the selection of the region of interest (ROI), Preset -allowing for the selection of a preset used for volume rendering and based on the type of the radiological data used (CT-AAA was used for CT scans of the nasal cavity), Shift -a slider button which determines the amount of voxels shown and which can be modified until the desired tissue is visible and as little of the undesired anatomy as possible is shown, Crop -enables the creation of a region of interest, Display ROI -visualises the selected region of interest, Rendering -VTK GPU Ray Casting should be used as the rendering option. After the volume is rendered and the region of interest selected, a new cropped volume should be created from the selected region. For this task, the Crop Volume Module is applied. The cropped volume and the region of interest in different views of the CT can be seen in Fig. 21.
The final step of the procedure performed in 3D Slicer involves the segmentation of the rendered volume. Segmentation is performed using the Segment Editor Module, in which the cropped volume from the previous step is used as the Master Volume of the Segment Editor. This module allows for the creation of a segment of the whole volume consisting of the tissue determined by a range of threshold values. These values should be used to determine the tissue to be displayed and the HU values for the desired tissue should be used as basis for the threshold range. As we are performing a CFD analysis, we want to encompass the volume of the nasal cavity which is filled with air, with as little of the tissue from which the internal walls of the nasal cavity are formed or the mucus, i.e. the fluid lining of internal nasal walls. In our example, we used the range from -1024 to -115, as this provided the best results. The resulting volume segment is shown in Fig. 22. Two more important properties of the Segment Editor Module are the Islands and the Smoothing menus. The Keep largest island options in the Islands menu removes any clusters of voxels (islands) which are not connected to the largest volume. The Smoothing menu can be used for producing a smoother surface of the volume segment. For the segmentation of the nasal cavity Median was used as the smoothing method with a Kernel size of 1 mm, but exact setting depend on the anatomy of interest. Lastly, the Segmentations Module should be used, with the Export and Models options selected, to export the selected segment as a separate volume segment. The Save Data icon can be used to save the volume segment as an STL file.

B.2. Model Clean-up and Preparation.
Once an STL model of the nasal cavity is obtained from the volume rendering and segmentation procedure, any freely available surface editing of 3D modelling software can be used to perform clean-up and final preparation of the model for the meshing procedure. The authors of this paper decided to use Blender [19] for this procedure, as it is freely available and has all of the features needed for 3D modelling. As there are many different applications available which are used for 3D modelling and most CFD engineers are comfortable in using at least one of them, the authors will not go into the specifics of editing and cleaning the model in Blender. Instead, an overview of the steps needed, independent of the application used, will be given next. Figure 23 depicts the model as exported from 3D Slicer. As the area in front of the nose is surrounded by air, it was also included as part of the STL model. What is more, there may be some leftover clusters of unconnected tissue that need to be removed. These issues need to be addressed during the clean-up phase of modelling, which includes. Once the box in front of the nose is removed and the nostrils exposed, the model needs to be separated into inlet and outlet patches. For the inlet domain, a spherical patch should be placed around the nose, encompassing the nostrils and allowing for later use of boundary conditions for mimicking atmospheric conditions. Finally, an outlet patch should be created near the exit of the nasal cavity in the direction of the trachea. The outlet patch should be created with a normal parallel to the direction of the main axis denoting the "down" direction to allow for easier application of the appropriate boundary condition. The cleaned model with separated inlet/outlet patches and prepared for the meshing procedure is shown in Fig. 23. B.3. Meshing Procedure. A step-by-step description of the meshing procedure will not be provided as the entirety of the meshing process is performed using cfMesh [16] -an open-source library for mesh generation, shipped as part of foam-extend, with tutorials and a user manual available online.
The model shown in Fig. 24 is used as the input model for the meshing procedure. The final model is comprised of three distinct patches shown in Fig. 2: cavity, trachea and atmosphere. Each patch is provided as a separate surface object in the STL file, simplifying the meshing procedure. The cartesianMesh algorithm is used, resulting in a unstructured mesh shown in Fig. 3. The cartesianMesh algorithm is based on inside-out meshing procedure, which requires a triangulated surface (e.g. STL file) to define the domain and to be used as a mesh template. The mesh is generated automatically based on the settings in the meshDict file and the provided input surface. The aforementioned meshDict file is provided in the case files distributed as part of this paper.
The generated computational domain is made out of predominantly hexahedral cells, resulting in a mesh with a nearly structured appearance. What is more, polyhedral cells are used as transition layers between regions of differently sized cells.   The simplified geometry was inspired by the geometry used in [2] and was created as a rough representation of the nasal passages, including turbinates. A representation of the simplified geometry, including boundary patches, is shown in Fig. 25. The names of the boundary patches were chosen as to be identical to the patch names of the real geometry used in this study. The computational domain was created using blockMesh, which is a basic mesh generation tool available as part of foam-extend. The resulting mesh is a block-structured mesh comprised of 11840 hexahedra. The computational domain based on the simplified geometry is shown in Fig. 26, while a cross-section of the right-hand side nasal passage is shown in Fig. 27.