Tomo-PIV in a patient-specific model of human nasal cavities: a methodological approach

The human nose serves as the primary gateway for air entering the respiratory system and plays a vital role in breathing. Nasal breathing difficulties are a significant health concern, leading to substantial healthcare costs for patients. Understanding nasal airflow dynamics is crucial for comprehending respiratory mechanisms. This article presents a detailed study using tomo-Particle Image Velocimetry (PIV) to investigate nasal airflow dynamics while addressing its accuracy. Embedded in the OpenNose project, the work described aims to provide a validation basis for different numerical approaches to upper airway flow. The study includes the manufacturing of a transparent silicone model based on a clinical CT scan, refractive index matching to minimize optical distortions, and precise flow rate adjustments based on physiological breathing cycles. This method allows for spatial high-resolution investigations in different regions of interest within the nasopharynx during various phases of the breathing cycle. The results demonstrate the accuracy of the investigations, enabling detailed analysis of flow structures and gradients. This spatial high-resolution tomo-PIV approach provides valuable insights into the complex flow phenomena occurring during the physiological breathing cycle in the nasopharynx. The study’s findings contribute to advancements in non-free-of-sight experimental flow investigation of complex cavities under nearly realistic conditions. Furthermore, reliable and accurate experimental data is crucial for properly validating numerical approaches that compute this patient-specific flow for clinical purposes.


Introduction
The human nose is a remarkable organ that serves as the primary gateway for air to enter our respiratory system.It plays a pivotal role in breathing, which involves oxygen intake and carbon dioxide exhalation.Further, the upper respiratory tract has important functions such as moisturizing, warming, filtering inhaled air, and detecting smell [1].Thus, pathological changes in this complex geometry and their impact on respiratory homeostasis significantly impair individuals' general wellbeing.Nasal breathing difficulties represent a severe medical concern, and their treatment is one of the most frequent surgical procedures performed by otorhinolaryngologists.In the United States alone, over ten million individuals are affected annually.This leads to substantial financial burdens on the healthcare system, causing an estimated annual cost of 22 billion dollars [2].
To reduce this substantial economic burden comprehensive strategies focusing on prevention, early detection, and costeffective treatment options are necessary.As such, understanding the dynamics of airflow within nasal passages is crucial for comprehending the intricate mechanisms involved in respiration, as well as various physiological and pathological conditions related to nasal function.In addition, medical procedures require a solid database for the planning of interventions.
The measurement of airflow within the nasal passages encompasses a broad spectrum of scientific research and clinical investigations.Numerous disciplines, including physiology [3][4][5] and biomedical engineering [6][7][8], have contributed to the advancement of our understanding of nasal flow dynamics.Flow measurements in the human nose aim to quantify the rate and distribution of airflow during inspiration and expiration, shedding light on factors such as nasal resistance [9][10][11] and the role of nasal geometry in shaping airflow patterns [12].Nasal airway geometry [13] and mucus viscosity [14] can influence the ease and efficiency of breathing.Techniques such as rhinomanometry directly measure the nasal airflow resistance by assessing pressure differentials across the nasal cavity [15].These methods provide information for diagnosing and monitoring nasal disorders, such as nasal congestion, septal deviations, and nasal valve collapse.
However, these methods fail to explain paradoxical phenomena such as empty nose syndrome (ENS) pathology.Furthermore, conventional flow measurement methods in the nose are insufficient for use as part of pre-operative planning to simulate the potential treatment outcome.
Understanding paradoxical flow phenomena in the nose requires more advanced imaging techniques, that integrate anatomical data from computer tomographic (CT) scans, such as Computational Fluid Dynamics (CFD) and Particle Image Velocimetry (PIV) to calculate and visualize airflow patterns within the nasal passages [8,12,16,17].CFD offers a complementary approach for studying nasal airflow by mathematical modeling and numerical simulation of fluid flow and transport phenomena in the nasal airways.Nevertheless, the suitability of these CFD investigations for implementation in the medical field requires careful consideration.Possible problems range from simplified geometries, neglecting the effects of cyclic respiratory flow, to the application of a suitable computational model [18].Further, different numerical approaches led to different solutions.As part of the OpenNose project, the approach described in this paper attempts to provide a basis for validation.
In contrast, the volumetric experimental investigation of PIV, tomographic PIV (tomo-PIV), provides valuable insights into complex three-dimensional flow phenomena, such as unsteady, bypass, and vortical flows, occurring during inspiration and expiration.By manipulating phantom models concerning nasal geometries and airflow conditions, researchers can study the effects of anatomical variations, nasal obstructions, and surgical interventions on nasal airflow patterns [6,19].Recent studies have already investigated the flows in the respiratory passages [20].The works range from stereoscopic (stereo) PIV investigations in the nasal cavity [6,21] with constant flow rates and flow rates based on breathing cycles [22].Researchers have utilized PIV to investigate alterations in airflow patterns caused by nasal disorders, such as septal deviations and high flow therapy [21][22][23].By providing PIV's boundary conditions and their results to CFD approaches, researchers can validate and refine computational models to enhance the accuracy and reliability of simulated nasal airflow patterns [6,[24][25][26][27].Recent numerical results are validated through simplified methods like fixed reference pressures in the nasal geometry.Further validation approaches include optical flow measurement techniques like stereo-PIV and tomo-PIV with simplified geometries or neglecting physiological boundary conditions such as the respiratory cycle.This integrated approach helps to bridge the gap between numerical predictions and experimental observations, providing a more holistic understanding of nasal airflow dynamics.
While PIV is a valuable tool for investigating nasal airflow and for creating a validation base for numerical investigation, a non-free-of-sight PIV investigation under realistic conditions poses unique challenges.One major challenge is the lack of physiological conditions: previous PIV experiments typically involve steady-state or simplified flow conditions to facilitate data acquisition and analysis.However, nasal airflow is influenced by various physiological factors, such as breathing patterns, anatomical structure, temperature, and mucosal properties.These factors can significantly affect the flow behavior in the nasal cavity, and their exclusion in PIV experiments may limit the realism and applicability of the findings to physiological conditions.Another challenge is the temporal and spatial resolution of the measurements: PIV measurements provide instantaneous velocity fields, capturing the flow patterns at a specific moment in time.Nevertheless, nasal airflow is highly dynamic based on the physiological breathing cycle and can exhibit significant temporal variations.The experimental setups of previous PIV measurements may not adequately capture these temporal variations, limiting the understanding of time-dependent flow phenomena.Additionally, the spatial resolution of PIV measurements may not capture fine-scale flow features, particularly in regions with complex geometries or small-scale structures within the nasal cavity.In addition, PIV experiments require clear and unobstructed optical access to flow fields.However, the nasal cavity's complex geometry, the presence of turbinate structures, and irregular surfaces can pose challenges in achieving undisturbed optical access owing to refractive index (RI) differences between the phantom models' material and the surrounding media.RI matching is a crucial technique to overcome these limitations by precisely matching the surrounding fluid's RI to that of solid objects (phantom models).Nonetheless, adjusting the surrounding medium regarding the RI changes the density ratio between PIV particles and the medium.This leads to an additional challenge: The PIV method assumes that the PIV particles are density-matched concerning the medium and follow the flow.
These problems are overcome by the hereinafter described tomo-PIV setup and investigations.This study provides highspatial resolved flow data in several regions of interest (RoIs) of a test person's nasopharynx considering different breathing cycle states.The combination with the micro-CT scan of the geometry will provide an experimental validation base for numerical investigations.The highly complex geometry combined with a cyclic flow requires the investigation of PIV parameters, such as the cameras' inter-double frame rate and the seeding strategies.Furthermore, the precise time control of the setup's components combined with repeated measurements provides analyses of unsteady flow phenomena via phaselocked investigations.Additionally, the study provides information regarding the reliability and accuracy of the experimental analysis.This project addresses methods to quantify accuracy such as the choice of PIV particles and the impact of ghost particles.Parameters such as the signal-to-noise ratio (SNR), particles' relaxation time, and standard deviations of repeated measurements are conducted.This procedure allows for a less biased patient-specific analysis of the nasopharynx.The described methodology allows the application of the technique to highly complex flows in human cavities.

Methods
Optical visualizations of flows within complex geometries necessitate dedicated experimental preparation, execution, post-processing, and finally a quantification of the accuracy and reproducibility of the experiments.After conducting the experiment, the methods for quantifying the accuracy of the experiment are discussed.
The following flow chart (figure 1) illustrates the sequence of necessary experimental steps.Figure 1 gives a summary of the necessary experimental steps, which will be described in the following paragraphs.The experimental procedure can be separated into the initial work (figure 1, red boxes), the preparatory work for the PIV experiment (figure 1, green boxes), the main PIV experiment (figure 1, orange boxes), and the post-processing work (figure 1, yellow boxes).The initial work (figure 1, red boxes) of the main process steps (figure 1, left) includes the manufacturing of the transparent silicone model followed by the RI matching (RIM) procedure between the silicone model and the surrounding fluid and the initial plate calibration of the tomo-PIV setup.The following preparatory work for the PIV experiment (figure 1, green boxes) considers the limitations of the double-frame camera system.To overcome the limitation of the small cameras' chip size, the large geometry of the transparent silicone model must be subdivided into smaller RoIs.The complexity of the nasopharynx coupled with the strongly varying induced flow rates based on the physiological breathing cycle, leads to another challenge: significantly diverse velocity magnitudes depending on the RoI and the phase of the breathing cycle (breathing cycle phase).As such, the cameras' inter-doubleframe time (PIV-dt) necessitates special attention.In general, Tomo-PIV conducts flow via particle motion between two successive images.Therefore, the particles' displacement within two frames (particle shift) strongly determines the accuracy of the reconstruction.As such, a single PIV-dt would not be sufficient to cover the entire breathing cycle.To overcome this issue, the particles' shift must be adjusted by the PIV-dt.In general, the recording frequency (limited to 15 Hz by the illuminating laser) is kept constant throughout the entire cycle.The breathing cycle is divided into equally spaced time slots and normalized by the entire length of the cycle.These are called breathing cycle phases (BCP).The PIV-dt was set for each BCP individually, to achieve a particle displacement of 4-8 pix in order to obtain reproducible measurements.The maximum allowable displacement was chosen on the basis of [28], where a maximum particle displacement of 10 pix is targeted for PIV measurement, and moreover, based on the postprocessing software providing a pre-defined setting for calculating a maximum expected displacement of 8 pix.The minimum value was set to 4 pix to ensure that accuracy is not compromised, according to the mean particle diameter of 3 pix.In other words, the PIV-dt has been adjusted according to the highly varying velocity fields due to the different flow rates based on a BCP to achieve a constant particle shift of 4 to 8 pix for the entire cycle.The seeding strategy is adjusted accordingly.After the preparatory work (figure 1, green boxes), the main PIV experiment (figure 1, orange boxes), including 100 repeated tomo-PIV measurements to establish a statistically-based data foundation, can be started.The RI must be controlled between repeated PIV measurements.The final steps include image pre-processing, tomographic reconstruction, and vector calculation (figure 1, yellow boxes).The interpretation of results is based on the arithmetic mean and the associated standard deviations.

Experimental procedure
Figure 2 shows a photograph (figure 2, left) and a schematic of the experimental setup (figure 2, right).In this context, the transparent silicone model constitutes the central basis of the experiment.Building upon it, the setup can be divided into three parts: The tomo-PIV part, the RIM part, and the flowinducing part.
2.1.1.Transparent silicone model.The creation of an optically transparent model, whose surface roughness and deviation from the physiological original can be considered negligibly small, requires the following steps: printing the inner geometry (3D rapid prototype), pouring the mold with silicone, and finally dissolving the inner geometry (figure 3).The basis for the transparent silicone model is a clinical CT scan of a test person's head (ethical approved) (figure 3(a)).To investigate the smallest cavities of the upper airways in a sufficient resolution an enlargement of the geometry is necessary.The post-processed, twice enlarged scan (magnification factor = 2), is then prepared for rapid prototyping (figure 3(b)).Given the large scale of the print, the trachea, as well as parts of the ethmoidal sinuses, are printed separately and fixed to the main print-body during coating (figure 3(c)).The prototyped model is made from gypsum powder and a binder (zp 131 powder/ zb 60 binder, Z Corporation, Burlington, USA) printed by ZPrinter 660Pro (3D Systems, South Carolina, USA).To stabilize it, rods with a diameter of 3 mm are placed at the sinus, at locations outside the line of sight of the measurements.Additionally, the rods minimized the final manufacturing step of the transparent silicone model (dissolution of the rapid prototype) by enabling better access to the inner geometry for the flushing water.The additional accesses allow for a higher seeding density in the subsequent PIV experiments and avoid separation of the working fluids component in the sparse flow-through areas as well.
The prototype acts as a negative for the silicone geometry (Sylgard 184, Dow Corning, Michigan, USA).The surface has to meet two criteria: it first, needs a minimum acceptable roughness and, second, requires a solid surface barrier to prevent the silicone from permeating into the negative rapid-prototype model.The two-component rapidprototype model is coated with 20 layers of a 2.5 wt vol −1 % polyvinylalcohol-distilled water solution (VWR chemicals, Pennsylvania, USA).The coating-sheet's thickness of the gypsum geometry is analyzed by a scanning electron microscope (Olympus Life Science Solutions, Hamburg, Germany).The processed rapid prototype is then placed in a rectangular  2.1.2.Tomo-PIV setup.The centrally placed fish tank, filled with the working fluid, builds the tomo-PIV setup base.It comprises a tomo-PIV system (LaVision GmbH, Göttingen, Germany), which includes three cameras with CCD sensors (Imager Pro X 2M, 1600 × 1200 pixels, LaVision GmbH, Göttingen, Germany).The cameras are equipped with macro lenses (ZEISS, Milvus 2/100), Scheimpflug adapters, and cut-off filters for better sensor protection.A double pulsed Nd:YAG laser (EverGreen200 70-200 mJ @ 532 nm, 15 Hz) is utilized.The laser is arranged perpendicular to the cameras' lines of sight (y-direction in figure 2).Knife edges are used to cut the oval-flared laser into a rectangular volume, the RoI.A RoI is defined by the initial scale of the field of view and the depth of field (L x × L y × L z = 49 × 36 × 6 mm 3 ) and is limited by the camera position setting.The flow is represented by fluorescence tracer particles to avoid, combined with cameras' cut-off filters, reflections of the surrounding light.The tracer particles (PMMA particles with Rhodamine B coating, with absorption resp.emission wavelengths of 560 resp.584 nm, LaVision, Goettingen, Germany) have diameters of 20-50 µm and a density ρ f of 1.19 g cm −3 .The tracer particles are induced either at the trachea or at the nostrils via two syringe pumps.All setup's components are centrally controlled using the software DaVis 10.2 (LaVision GmbH, Göttingen, Germany).The phantom model is then placed in an hexagonshaped fish tank containing glass plates, oriented at a 60 • angle to each other.The alignment of the glass plates allows for an orientation of the three cameras' lines of sight at a 60 • angle, perpendicular to the particular glass plate.This minimizes the disadvantage of the optical distortions inherent to the reflections caused by the glass plates' surface.A 60 • angle is chosen to ensure both, static stability of the fish tanks' construction and a satisfactory tomo-PIV reconstruction quality [29].

Image capturing and preprocessing.
Images are captured by three CCD cameras.The capturing mode is double-frame.To achieve an average particle image shift of 4 to 8 pix between two consecutive images, the PIV-dt is preliminarily set in a range of 2000 µs-15 000 µs, based on the breathing cycle phase of interest.The captured images were pre-processed to achieve a background pixel intensity of 0 cts.Furthermore, since the complex geometry includes nonseeded areas (only silicone) and seeded regions, a mask is applied to each viewing angle in the x-y slide.For calibration purposes, the first mapping function of the three cameras' location in a 3D space is calculated using a plate calibration target.The calibration is further improved by volume self-calibration (VSC) [30].The size of the image reconstructed by the mapping function is 52 × 40 mm 2 (final RoI L x × L y × L z = 49 × 36 × 6 mm 3 ) with a scale factor of 29 pix mm −1 .For the VSC images, an average seeding density of 0.006 particles per pixel (ppp) in the same experimental setup is used.

Reconstruction and vector calculation.
The image reconstruction and the vector calculation are prepared using the software DaVis 10.2.Following to image pre-processing, a volumetric reconstruction is fulfilled by a first initialization using the MLOS algorithm, followed by a MART algorithm with six iterations.Then, sparsity determination is fulfilled, followed by two iterations of the SMART algorithm.The vector calculation is processed for the interrogation window with a maximum expected displacement of eight voxels and a final window size of 32 voxels.

RI matching.
The two components of the working fluid are of different densities.This leads to a segregation of the working fluid components, particularly in the 'low flow regimes' of the geometry.Thus, an additional effect occurs besides evaporation, which causes the RI to shift.Online monitoring of the RI is necessary to detect this shift at an early stage and to intervene between two measurement recordings.As such, a new method of high-accuracy RIM is developed.To achieve a match between the RI of the phantom model and the working fluid (water-glycerol-solution), the mixing ratio between the two components of the working fluid is adjusted using the cross-line laser setup shown in figure 2. Here, the red planes of the laser-light-cross spanned by cylindrical lenses traverse the measurement volume horizontally.The laser-light cross crosses the investigated RoI and is visualized on a scaled paper upon leaving the fish tank.Images of this visualization are captured by a camera equipped with a CMOS chip of 6000 × 4000 pix (D3400, Nikon Corporation, Tokyo, Japan), and a magnifying lens (AF-P DX NIKKOR, 22/18, Nikon Corporation, Tokyo, Japan).First, a reference image of the laser-light-cross without the model is captured.After placing the phantom model in the fish tank further images are captured and the RI of the working fluid is measured using a digital refractometer (ORM 1RS, KERN & SOHN GmbH, Balingen, Germany).In order to adjust the RI of the working fluid, and finally eliminate optical distortions leading to noisy tomo-PIV results, specific amounts of the working fluid components are added and mixed.Meanwhile, the movement of the laser-light cross is monitored and the RI is measured after each change in the mixing ratio.

Breathing-cycle-based flow rates.
The pump is controlled using a programmable logic controller (PLC) and TwinCAT software (Beckhoff Automation GmbH, Verl, Germany).To investigate the flow field, 30 breathing cycles from a healthy subject are recorded using a spirometer (Pneumotrac TM , Vitalograph @ GmbH, Hamburg, Germany), and processed with proprietary Spirotrac 6 software.The flow rate data are adjusted using dynamic similarity to match the experimental conditions.Therefore, Reynolds' similarity is ensured based on equation (1) with a magnification factor M, the experimental flow rates Q vitro , the spirometric flow rates Q vivo , and the corresponding dynamic viscosities η vitro and η vivo .
The dynamic viscosity of the working fluid η vitro is determined using a rotational rheometer (Discovery HR20, TA Instruments, New Castle, USA).MATLAB software is used to process the flow rate cycles with dynamic time-warping and fit them with a 7th degree polynomial function.Womersley similarity is ensured based on equation (2).Therefore, the overall length of the experimental breathing cycle ∆t vitro is calculated from the dynamic viscosities η vitro and η vivo and the overall length of the experimentally captured cycles ∆t vivo .
2.1.5.Device control and phase averaging.All setup components are centrally controlled by a master trigger, which induces the flow rates by initializing the linear motor-driven piston pump (figure 4).For lowering experimental costs, a seeding strategy was necessary.Thereby, two main challenges arose: On inspiration, particles added to the trachea would only enter the pump.On the other hand, during exhalation, the particles that were added at the nostrils would be expelled into the environment.The seeding strategy was set up based on the RoI and the breathing cycle phase.As such, based on the master trigger, the two seeding pumps are controlled, either inducing particles in the trachea or in the nostrils.Therefore, the seeding pumps' speed, their initial seeding time, and their total seeding time are adjusted preliminary regarding the investigated breathing cycle phase and the RoI individually for each pump.Following the precise temporal device control, phaseaveraging methods enable the identification of recurring complex flow phenomena as well as non-steady, chaotic flow phenomena.In addition to the accurate and reproducible representation of the flow situation during a defined phase of the breathing cycle, depicting various BCPs poses an additional challenge.Therefore, especially cameras with CCD sensors are limiting.In addition to the general advantage of CCD sensors in terms of improved SNR, CCD chips have the drawback of being 'low-speed' systems.In low-speed systems, a single PIV-dt must be predetermined and remain constant throughout the entire recording.Consequently, the initial selection of the breathing cycle phase of interest becomes dependent.Certain phases of the breathing cycle require different PIV-dts.Hence, only a few cycle phases can be covered within a single measurement.Through repeated measurements and the adjustment of the PIV-dt, various cycle phases can be analyzed for several RoIs.To achieve this, results are analyzed within 100 repeating cycles.A breathing cycle phase is a fixed time frame within a breathing cycle.Different BCPs can be accessed by adapting PIV-dt from 2 ms up to 15 ms.

Accuracy investigations
For the 3D reconstruction of tomo-PIV, the a-priori setup parts need to be very accurate [28].The following paragraphs describe methods to directly address the particles' ability to follow the flow as well as the accuracy of the reconstruction of the experiment.

Particles and working fluid.
The tomo-PIV method relies on illuminating and visualizing tracer particles added to the flow.This requires a negligible influence of inertial and gravitational forces on the particles.As such, the error introduced by the tracer particles' density and gravitational forces, as well as the tracer particles' ability to follow the flow [31] is investigated by assessing a density ratio, the Stokes law and by their relaxation time [28,31,32].
Gravitational forces can cause errors, if the density of the working fluid ρ f and the particles' density ρ p are unequal.Therefore, it is essential to ensure a nearly equal density between the working fluid and particles.The density of the working fluid (RI = 1.4139) is approximately ρ f ≈ 1.15 g cm −3 .With a particle's density ρ p of 1.19 g cm −3 , the ratio between the densities ρ fp of approximately 0.98 is calculated based on Oucheggou et al [33].
Investigations on tracer particles ability to follow the flow are carried out on an exemplary RoI (at the trachea): The hydraulic diameter d h is 25 mm in this exemplary measurement region, and the average velocity magnitude of the fluid is set to v ≈ 0.08 m s −1 .The estimations are made for the largest particles with a diameter of d p = 50 µm since smaller particles will lower the relaxation time [34,35].The Stokes law is valid for a creeping flow around a single particle and small fluid's Reynolds numbers Re f .As such, the Reynolds number for particles Re p < 1 (diameter d p > 1 µm), Re f , and the necessary Stokes velocity of the particles v St have to be addressed for further investigations [31][32][33].The working fluid's dynamic viscosity η f is, due to RIM and the amount of the two components, 9.44 mPa s. v St (particle diameter d p set to 50 µm) is 1.153 • 10 −5 m s −1 .Re p is 7.02 • 10 −5 and the estimated Re f is 2436.Now, assuming the Stokes law to be valid for small fluid's Reynolds numbers Re f and Re p <<1, the Stokes number is approximately 0.03 and further, the particles' relaxation time τ p is approximately 1.7 • 10 −6 .

Reconstruction accuracy.
To obtain an overview of the quality of the reconstruction and calculation for such a complex geometry, the light intensity signal-to-noise ratio SNR E and the particle distribution signal-to-noise ratio SNR N for all intensities as well as for the real particle intensities SNR PI170 are investigated.Investigations are conducted on the reconstructed intensities of the calculated volume (L x × L y × L z = 51 × 51 × 60 mm 3 ) of an exemplary double frame.The sum of the reconstructed intensities of all y-slides is displayed in x-z-projection using DaVis software.Increasing the reconstruction and calculation field well outside the illuminated region enables the calculation of the ratio SNR E , which correlates the average intensity value inside the illuminated volume E I and that outside the illuminated volume E O [36].
Further, the particle distribution signal-to-noise ratio SNR N is investigated with the total counted number of actual tracer particles N P in the tomographic reconstruction and the total counted number of ghost particles N G .
The counted number of particles inside the illuminated volume per z-axis value (per camera's chip pixel row) results from the number of ghost particles plus the number of actual particles (real particles) [36][37][38].A histogram, representing the number of pixels per z-axis value and their intensity values inside and outside the light volume per z-axis value, can be transferred into the probability density function of the ghost particles' peak intensity values.The values inside the illuminated volume minus those outside the illuminated volume provide an estimate of the real particle peak intensity distribution.This enables determining the ratio of N P and N G at the average peak intensity of actual the particles SNR PI [39].

Results
The methodology for investigating complex internal flows in the upper human airways allows the visualization and interpretation of the intricate flow patterns in the nasopharynx.
Creating a transparent silicone model based on a clinical CT scan, precise matching of refractive indices, and timecontrolled measurements reveal complex and recurrent flow patterns that are dependent on the geometry during a physiological breathing cycle with quantifiable accuracy.

Transparent silicone model
The precise execution of the manufacturing steps while adhering to the defined parameters, including manufacturing a rapid prototype based on a clinical CT scan, surface processing of the rapid prototype, model casting, and dissolution of the rapid prototype, enables the production of a transparent and optically accessible silicone model.During the surface processing step, the non-coated sample of the prototype exhibits an average surface roughness of 9.707 µm at the first and 7.534 µm at the second measuring point.After coating, surface roughness values of 0.074 µm respectively and 0.099 µm are achieved.The sheet thickness reaches 0.18 mm compared to the overall wall thickness of the 3D-printed negative (3 mm).
Investigations on the surface roughness pre-and post-coating show an average decrease of nearly 98% of the surface roughness while only increasing the wall thickness by 6%.As a result, the surface quality of the 3D print can be transferred to that of the phantom model.

Image capturing and preprocessing
The initial plate calibration results in an error of 0.34 pix.By VSC with a final sub-volumes' number of 10:10:4 (x:y:z), an average error lower than 0.1 voxels is reached.Image pre-processing increased the signal-to-background ratio from 500 cts.(signal particles) vs. 200 cts.(signal background) to 1500 cts. vs. 0 cts.The processing results in particle diameters between 1.7 and 3 pix.

RI matching
The working fluid consists of two components with their RIs, water n(H 2 O) = 1.333 and glycerol n(Glyc) = 1.473.As shown in figure 5, the described method enables a sensitive RIM by analyzing the disruption as well as the displacement of the projected lines of the laser-light-cross on the scaled paper (figure 5, black cross, laser-light-cross not passing the  phantom head).Figure 5 shows the visualizations of the laserlight-cross while passing the phantom head and the working fluid comprising different RIs (1.4135-1.4145)(figure 5, color-coded crosses).A RI of 1.4139 results in solid lines and a laser-light-cross matching the reference.A RI of 1.4139 is valid for a mixture comprising 52 parts glycerol, leading to a VSC error lower than 0.1 pix.Shifts in the amount ratio of the components of the working fluid result in spreading crosslines and a shift of the cross compared to the reference image.Furthermore, online monitoring resulted in satisfactory VSC's over longer periods of measuring.

Breathing-cycle based flow rates
With the flow-inducing part, the flow rates are induced by a linear motor-driven piston pump at a certain flow rate, based on test person data adjusted by the Reynolds and the Womersley similarity.Figure 6 visualizes the flow rates' curve development induced by the pump.While the modification by Reynolds similarity based on equation ( 1) requires a factor M•(ν vitro /ν vivo ) = 0.98 for the flow rates, the modification by the Womersley similarity based on equation ( 2) requires a factor of M 2 •(ν vivo /ν vitro ) = 8.9.The final flow rates cover a cycle period of 24.8 s, including an expiration part of 13.72 s.Thereby, the peak flow rates are achieved at 6.0 s (0.370 l s −1 ) for expiration and 17.7 s (−0.372 l s −1 ) for inspiration.All components are centrally controlled using a PLC.The flow is introduced into the trachea and leaves the geometry at the nostrils during inspiration.

Phase locked processing
The following presents the flow patterns throughout the entire respiratory cycle for one RoI inside the nostril and one RoI at the end of the trachea where it transitions into the nasal geometry.Figure 7        Velocity magnitudes during the breathing curve's flow rate peaks reach up to 0.2 m s −1 .In this context, the exhalation phase exhibits higher absolute values of local velocity magnitudes.Strong differences in flow patterns between the inhalation and exhalation phases are highly noticeable.While the exhalation phase exhibits a circular flow structure formation before entering the conchae, the inhalation phase shows a highly directional flow.When examining the constant induced flow rates' flow patterns with the respiratory cycle-based flow rates' flow patterns, the pattern at BCP 0.98 (figure 9(n)) closely resembles that of the CFR (31.4 ml s −1 ) (figure 9(b)).The velocity magnitudes throughout the RoI are nearly identical.However, the vectors in the upper left corner of the geometry differ significantly, exhibiting an unstructured flow direction.Contrary to the RoI located at the nostril, the flow pattern located at the trachea during the constant exhalation flow rate (−31.4 ml s −1 ) (figure 9(a)) noticeably differs from the flow pattern at BCP 0.52 (figure 9(h)) in the respiratory cycle.While the result of cyclic volume flows exhibits higher velocity over most of the geometry, the constant measurement shows a more pronounced jet.
The precise temporal control of the measurements allows for examining the evolution of the standard deviations of phase-locked velocity profiles throughout the entire cycle.This enables drawing conclusions about recurrent flow phenomena over time as well as changes in flow phenomena within the same phase across repetitive cycles.Figure 10 shows the mean velocity (cross markers) across the entire RoI (Nostril in red, Trachea in blue) over the breathing cycle.Additionally, the standard deviation calculated from the arithmetic mean of a single breathing cycle phase for each RoI (over 100 repeating cycles) is shown throughout the breathing cycle (circle markers).The breathing cycle curve is clearly discernible throughout the course of the mean velocity magnitude.In this regard, the graphs representing the RoI located in the nostril and in the trachea exhibit similar trends.During the transition from exhalation to inhalation part, low mean velocity magnitudes are observed (<0.01 m s −1 ), while the transition from inhalation to exhalation part has mean velocities higher than 0.02 m s −1 .The transition from exhalation to inhalation part occurs, corresponding to the breathing cycle curve, between 13 and 14 s.The curve of the nostril's RoI exhibits a slight phase shift compared to that of the trachea at later time points.The plotted standard deviations of the phaselocked velocity averaging process show values ranging from 0.015 m s −1 to 0.025 m s −1 .The standard deviations exhibit increasing values with increasing mean velocity magnitudes.While the standard deviations at the trachea generally show higher values during the inhalation phases, the standard deviations at the nostril indicate higher values during the exhalation phases.

Reconstruction accuracy
Reconstruction accuracy investigations are conducted on the reconstructed intensities of the calculated volume of an exemplary double frame.The sum of the reconstructed intensities of all the y-slides is displayed in the x-z-projection.The calculation in a rectangular area of 1370 × 245 pix results in a SNR of 2.3 between the averaged intensity well inside and well outside the illuminated volume (figure 11

Discussion
Conducting high spatial resolution experiments using tomo-PIV offers a comprehensive and detailed understanding of nasal flow, contributing to advancements in medical research and engineering applications.The study provides a threedimensional view of the flow field within parts of the nasal cavity, allowing researchers to understand the complex flow patterns in great detail.This research project considers the necessary physiological conditions.To achieve a high level of accuracy, the study conducts a detailed analysis of the measurement's reliability influencing parameters.In order to investigate the breathing cycle's impact on the induced flow rates, the patterns of the natural respiratory cycle-based flow rates and the flow patterns of a CFR are compared.
The two-fold magnification of the geometry enables accurate analysis of smaller flow structures.While constructing the phantom head based on real patient CT data, assessments of individual flow structures possible without neglecting flow-determining parts, such as geometrical influences as well as the physiological breathing cycle.The roughness of the rapid-prototype model shows a huge improvement in the surface finish at 98% while only having a negligible role in the change in the wall thickness (6%).The final surface roughness of the geometry is less than 0.4% of the smallest tracer particle diameter.The smooth surface of the model enables flow investigations at the geometry's near-wall regions, while not causing any errors, especially in smaller cavities.
The new method of RIM minimizes the remaining geometrically induced optical distortions.Directly integrating the RIM setup into the experimental setting enables real-time RI monitoring, precise RI adjustment, and consequently highprecision RI verification in a highly complex geometry.The RIM method associated with a laser light cross enables an average disparity, after VSC, of less than 0.04 voxel (despite the need to have the line of sight pass through four different types of media, each with a different refractive indicex (air, glass, working fluid, silicone)).
Nonetheless, PIV requires a nearly equal density of tracer particles and the surrounding media.However, the density of the working fluid and its RIs are determined by the ratio of its two components.A change in the RI also results in a change in the working fluid's density.It is crucial to precisely assess the accuracy of particles following the flow of the working fluid with a perfectly matched RI.Investigations on the particles revealed a small difference between the density of the working fluid and that of the tracer particles.Despite such small differences, the suitability of the particles with respect to the working fluid must be further investigated to achieve sufficient measurements.ρ fp = 0.98 nearly equal to one, gives a good indication for the matching working fluid's and the tracer particles' density.The very small Re p ≈ 7.02 • 10 −5 permits the assumption of a creeping fluid's flow around a tracer particle as well as the calculation of the Stokes number.St ≈ 0.03 indicates suitable tracer particles for St → 0. Thus, τ p ≈ 1.7 • 10 −6 s remains very small compared to the shortest choice of the PIV-dt = 2 • 10 −3 s.The relaxation time might end up being just about 0.08% of the analyzed time.From this assumption, a suitable combination of tracer particles and the working fluid can be inferred for this experimental setup in regard to the choice of the PIV-dt.
Image pre-processing resulted in the intensity of the background of 0 cts.with an average diameter of the particles between 1.7 and 3 pix.The 3rd-order polynomial function of the initial plate calibration in an error of 0.38, which remains below the recommended set at <0.4 pix [30].The best result of VSC achieved with sparse seeded images (0.006 ppp) and the reconstruction error can be minimized to a value <0.03 pix.Note, that the presented setup also an adjustment of the seeding density.Given the recommended values SNR E ⩾ 2 [40] we can achieve a value SNR E = 2.3 which is a first good assumption for a sufficient reconstruction result.Furthermore, the calculated SNR N of 1.2 is higher than the minimum recommended value of SNR N > 1 [36][37][38].The SNR with regards to the average intensity of the real particles, reached a value of SNR PI170 = 2.3.Besides the typical number for the uncertainty of PIV displacement calculation set at 0.1 pix units [41], SNRs showed sufficient values for a highresolution proof of concept measurement.
The methods described investigate the flow based on a physiological breathing cycle on exemplary positions of the phantom head.The goal of applying the experimentally collected data in a clinical setting necessitates a precise design, execution, and evaluation of the experiment.Since the experimentally investigated breathing cycle is modified by dynamic similarity, complex physiological flow phenomena are easy to investigate.Thus, the study replicates the physiological airflow patterns and allows interpreting of them on a medical base.The adjustment through Reynolds' and Womersley's similarity allows for representing the actual airflow in the airways during physiological breathing by using a working fluid in a geometry up by a factor of two.The influence of inertia forces and geometry becomes significant through the factor securing Womersley similarity (factor 8.9).Another significant advantage is the ability to perform phase-averaged evaluations.
The phase-locked averaging process, has proven to be an effective method for capturing the magnitudes and flow directions within a defined RoI.The phase-locked averaging process helps to filter out noisy measurements, which may be caused by the experimental setup including a variety of devices.However, using CCD cameras to investigate the flow patterns occurring during the physiological breathing cycle restricts the number of cycle phases that can be covered in a single measurement.This limitation necessitates the careful selection of the breathing phase of interest and the adjustment of PIV-dt for accurate measurements.Adjusting the PIVdt concerning the RoI and the breathing cycle phase allows for a detailed examination of flow patterns during specific phases of the respiratory cycle.Notably, the visualization of flow patterns during exhalation and inhalation reveals distinct differences, including the formation of a circular flow structure during exhalation and highly directional flow during inhalation at the trachea as well as at the nostril.A highly accurate proof of concept method for analyzing flow structures using 3C3D vector fields in human cavities is indispensable.The importance of a three-dimensional, spatially high-resolution examination of the flow patterns in the upper airways is evident in figure 8.While a two-dimensional analysis of the vectors indicates a circular structure, a three-dimensional examination reveals a spiral flow structure.
In addition, phase-locked averaging identifies unsteady flow phenomena as indicated by the increasing standard deviations when high average velocity magnitudes (figure 10).The minimum values of the standard deviations for both investigated RoIs are approximately 0.0025 m s −1 .This indicates the measurement uncertainty of the The effects of minimal time delays in the phase-locked process, ghost particles, and measurement inaccuracies of extreme velocity differences within a RoI become apparent.Another problem is the flow near the geometry's boundaries.Owing to the adjustment of PIV-dt to the velocity of the main flow, the flow near walls exhibits a particle shift that falls below the resolution limit.This is also evident when looking at local standard deviations within a RoI.The higher relative local standard deviations are particularly noticeable in the areas close to the walls.However, the higher standard deviations at higher mean velocities suggest that different flow phenomena occur in repeating cycles.This is particularly evident during exhalation in the RoI at the nostril compared to inhalation in the RoI at the trachea.The application location of the flow can explain the visible phase shift between the trachea and nostril.The flow is introduced at the beginning of the trachea.Due to inertia forces, flow phenomena initially occur closer to the application site before eventually reaching the nostril.

Conclusion
A high-accuracy proof of concept method for analyzing flow structures by 3C3D vector fields in human cavities, based on the patient's individual geometry, is possible through the presented construction of a realistic in − vitro model and a high-resolution visualization of the flow structure.Investigating the complex flow patterns during breathing at two locations in the nasopharynx sheds light on the complex network of flow structure determining boundary conditions like the induced flow rates and the geometry.The different flow structures in the dependency of the breathing cycle phase illustrate that boundary conditions, like breathing cyclebased flow rates, are indispensable to achieve meaningful results.Further, small spiral-shaped structure that occurs, as visible in the present study, influences the flow conditions in the entire complex geometry.It is conceivable that surgical interventions at the depicted position, which contribute to airflow through the turbinates, may lead to long-term consequences such as inadequate airflow through the lower nasal turbinates if the interventions are not performed correctly.Additionally, the method can be applied to flow analyses in other highly complex geometries, e.g.blood flow applications based on patient-specific CT scans, and offers a means of validating CFD datasets for potential implementation in the clinical environment.

Figure 1 .
Figure 1.Flow chart of the experimental procedure.The experimental procedure can be separated into the initial work (red boxes), preparatory work for the PIV experiment (green boxes), the main PIV experiment (orange boxes), and the post-processing work (yellow boxes).RoI individual steps are marked in purple, and the breathing cycle phase individual steps are marked in gray.

Figure 2 .
Figure 2. Scheme (left) and photograph (right) of the measurement setup, visualizing the the centrally placed phantom model.Around the fish tank, the tomo-PIV setup, the flow-inducing part, and the RIM part are placed.Three CCD-cameras with optical equipment, triggered by a programmable timing unit are included.Their line of sight is geared to the measurement region, which is illuminated by a light volume emitted by an Nd:YAG laser.The light volume is cropped by knife edges.Fluorescent tracer particles are added to the flow by two syringe pumps.A linear motor-driven piston pump generates the flow rates.

Figure 3 .
Figure 3. Performed processing steps during phantom head manufacturing, (a) data captured by clinical CT imaging, (b) prepared geometry for rapid prototyping, (c) printed model during coating, (d) casted model in a rectangular box, (e) molding procedure by filling silicon into the box, (f) de-molding with visible prototype in the airways, (g) de-molded phantom head indicating characteristic anatomical parts.

Figure 4 .
Figure 4.The hierarchical trigger protocol of the time-dependent setups component.The master trigger with a 5 V TTL signal initializes the breathing cycle comprising the breathing cycle phases (BCP), the seeding by two syringe pumps, as well as the cameras' capturing mode.

Figure 5 .
Figure 5. Results of the RIM method.The laser-light-cross passing the working fluid (black cross) and the phantom head (color-coded crosses) is visualized on a scaled paper.The disruption and the shift of the laser-light-cross caused by mismatched RIs of the phantom head model and the working fluid are clearly visible.

Figure 6 .
Figure 6.Development of the physiological breathing cycle based induced flow rates.Experimentally captured flow rates during repeating breathing cycles are modified and adapted by the Reynold's and the Womersley similarity.
2 s) to give the opportunity for comparing the flow situations with the continuous flow rates' (CFR's) situation.When selecting the phases, care was taken to adequately represent the different flow situations.In addition, flow patterns in the nostril during a constant pump's flow rate of 31.4 ml s −1 respectively −31.4 ml s −1 are visualized in figures 7(a) and (b).The constant flow rates were chosen in order to address physiological value included in the cycle and moreover, match boundary conditions of numerical investigations.The velocity magnitudes are color-coded, and the indicated vectors represent the directions of the local velocities.The increasing and subsequently decreasing magnitudes of the velocities corresponding to the rising and falling flow rates due to physiological breathing are clearly visible.While the peak velocity magnitudes, corresponding to phases of the cycle with low induced flow rates (points near the intersection of the breathing curve

Figure 7 .
Figure 7. Phase averaged mean velocity fields at the nostril at a constant pump flow rate (a), (b) as well as during various BCPs (c)-(n)).

Figure 8 .
Figure 8. Flow structures of the left nostril.Streamlets visualize flow structures in the investigated RoI in the left nostrils, beside the nasal septum.The frontal view (a) and (c) shows the circulating structure in a color-coded velocity magnitude.The rotated view (b) and (d) shows a helical structure of the streamlets.
Figure 9 illustrates the flow patterns occurring in the RoI located at the end of the trachea, entering the main nasal geometry.The same BCPs as those in figure 7 are depicted.The velocity magnitudes are color-coded, and the indicated vectors represent the directions of the local velocities.Furthermore, the occurring flow patterns during a CFR (31.4 ml s −1 , −31.4 ml s −1 ) (figures 9(a) and (b)) induced by the pump are illustrated.

Figure 9 .
Figure 9. Phase averaged mean velocity fields at the trachea at a constant pump flow rate (a), (b) as well as during various BCPs (c)-(n)).

Figure 10 .
Figure 10.Average velocity magnitude vs. standard deviations of the phase-locked average method.
left green box).The left part of figure 11 shows the total number of counted intensities per z-axis value (z-axis value corresponds to a pixel value of the camera chip) (figure 11 blue graph).Means for N G ≈ 7.5 of particles well outside the light sheet and N G + N P ≈ 16.5 particles per z-axis value (directly connected intensity values represent a particle) inside the illuminated volume are indicated in the probability functions in the right part of figure 11.These values result in an SNR N of 1.2.The SNR at the average ghost particle peak intensity, based on the probability density functions of the peak intensities, is calculated at SNR PI170 ≈ 2.3 as shown in the right part of figure 11.

Figure 11 .
Figure 11.Visualization of the accuracy investigations.The left part of the figure showcases the sum of y-slices'(of an exemplary double frame located at the trachea) intensities displayed in the x-z projection.The blue graph shows the number of counted particles in the z-axis, where the green lines indicate the illuminated part.The right part visualizes the calculated probability density function (pdf) of the different counted intensities.