Computational fluid dynamic models as tools to predict aerosol distribution in tracheobronchial airways

Aerosol and pollutants, in form of particulates 5–8 μm in main size face every day our respiratory system as natural suspension in air or forced to be inhaled as a coadjutant in a medical therapy for respiratory diseases. This inhalation happens in children to elderly, women and men, healthy or sick and disable people. In this paper we analyzed the inhalation of aerosol in conditions assimilable to the thermal therapy. We use a computational fluid dynamic 3D model to compute and visualize the trajectories of aerosol (3–7–10–25 µm) down to the sixth generation of bronchi, in a steady and dynamic condition (7 µm) set as breath cycle at rest. Results, compared to a set of milestone experimental studies published in literature, allow the comprehension of particles behavior during the inhalation from mouth to bronchi sixth generation, the visualization of jet at larynx constriction and vortices, in an averaged characteristic rigorous geometrical model including tracheal rings. Results on trajectories and deposition show the importance of the including transient physiological breath cycle on aerosol deposition analyses. Numerical and graphical results, may enable the design of medical devices and protocols to make the inhalations more effective in all the users’ population.


Scientific Reports
| (2021) 11:1109 | https://doi.org/10.1038/s41598-020-80241-0 www.nature.com/scientificreports/ compared data from simplified and realistic models of TB tree down to the 3rd generation with the empirical data of Zhou et al. 28 to quantify the influence of geometrical features (i.e. tracheal curvature, main bronchi curvature and irregular cross sections) on pharmaceutical aerosol deposition. They found that accurate tracheal features enhanced regional deposition and shifted particle patterns towards left walls. Realistic CSAs and bifurcation curvature increased DE in the TB tree 3 . Moreover respiratory air flow characterized in the lung is nonsymmetrical while the pattern of flow field was similar. Das et al. 20 used an idealized, anatomically-faithful upper airway geometry, where simulations are function of age, from a 5 year old to an adult. The results of the comparison between a Dry Powder Inhalers (DPI) and nebulizer inhalation performance were collected in a dimensionless curve governing deposition in the airways via Stokes number. A patient-specific model was then simulated by Farghadan et al. 10 within a breathing cycle, and the use of Finite-time Lyapunov exponent (FTLE) was involved to compute the Lagrangian topological maps that define the destination of particles at the trachea. Tian et al. 16 compared steady state and transition simulations of particles deposition in a characteristic patient model, by means of a Stochastic Individual Path (SIP) model in branches deeper than TB4.
CFD results were compared for validation to experimental data on total and regional DE. Experimental data were obtained from both human volunteers and airway hollow casts [28][29][30][31][34][35][36][37] . For example Chan and Lippmann 31 studied the regional DE in a hollow cast of the human larynx and TB tree extending down to the 6th generation, and in 26 human volunteers, finding a linear dependence of particle DE on Stokes number for aerosols with aerodynamic diameters greater than 2 µm.
Similarly, Cheng et al. 29,34 , investigated the effects of particle size and breathing conditions on DE in a human oral airway replica. They observed that DE increased with higher flow rates and particle diameters, suggesting that impaction is the dominant deposition mechanism. Furthermore, by comparing their results with in-vivo deposition data 31,37 , they defined a relationship between DE and the impaction parameter IP, defined as: where ρ is the density of the particle, d p is the aerodynamic diameter and Q is the flow rate 29 . Regional DE in a TB tree replica was carefully studied by Zhou et al. 28 through the injection of fluorescent particles at constant flow rates (15,30, 60 L/min). They observed a relationship between DE and Stokes numbers, bifurcation angles and diameters and provided an empirical model to estimate DE in each generation.
Many partial/complete model of human airways were create to investigate airflow pattern and particle DE under various breathing conditions. Airflow pattern was widely analysed under both constant flow rates 5,9,13,15,26 , and transient breathing cycle 14,18,38 , while most of particles analyses found in the literature were performed under constant airflow conditions. Gurman et al. 36 suggest that this approximation can represent large particle cases pretty well while it underestimates total DE of smaller particles.
Only a few recent works analysed the behaviour of nanometric and micrometric particles during a complete breath cycle 10,16,24,38 , under different breathing conditions (sedentary, light and heavy breathing).
The focus of the present research is to create a characteristic model of human airways, in order to investigate therapeutic aerosol DE, including the transient behaviour of a physiologically breathing cycle in sedentary conditions. A deep analysis of the anatomical features described in literature suggested the involvement of a special strategy to design the geometry of bronchi around bifurcations and the introduction of the cartilaginous rings in the reconstruction of the trachea. Presence of rings, usually considered as structural components, modifies the appearance of the tracheal wall as sub-millimetric features 16 that may have a role in the formation of boundary layer of flow and in the deposition of particles, consistently with in-vitro observations of Russo et al. 39 . The model of human airways has been constructed to allow future scalability based on the age-specific anatomy.

Results
Results of numerical simulations are tools to visualize and describe the behaviour of air and aerosol along a respiratory cycle.
Airflow. Ventilation was quantified considering the airflow exiting the RUL (right-upper lobe), the RLL (right-lower lobe) and the RML (right-medium lobe), while left ventilation took into account the LUL (leftupper lobe) and the LLL (left-lower lobe) outflows. Reasonable consistency of the percentage ventilation with published data [40][41][42] , was found for each lung lobes: RUL 19%, RLL 29%, RML 9%, LUL 26% and LLL 17%. The total left and right total ventilation values were 43% and 57%, respectively. The resultant distribution is translated to the transient simulation, by means of dedicated script. The velocity profile at each cross section varies continuously during one complete cycle of respiration, depending on cross sectional areas, the geometrical features and the velocity profile at the inlet of the region (Fig. 1a) 1 . As expected, a laryngeal jet appears in the glottis region, where the noticeable constriction of the larynx forces the airflow to accelerate 16,20 . The strong centrifugal forces shift this high velocity region towards the posterior wall during the inhalation phase (0-2 s), while the peak is shifted towards the anterior walls during the exhalation phase (2-4 s). Nevertheless, airflow always exhibits the characteristics of laminar or quasi-transitional flow, reaching a maximum Reynolds number equal to Re = 1848 at t = 1s in the glottis, which is the most critical section in terms of area and local velocity. The absence of a turbulent regime is reasonable considering the sedentary breathing condition: turbulent phenomena appear for airflow greater than 10-12 L/min 1,9,14,26 .
The static pressure distribution in Fig. 1b shows a pressure drop equal to 4.5 Pa between the tracheal inlet and the 6th branches at the inspiratory peak, which is close to measurement by Gemci et al. 9 .
The laryngeal jet effect propagates downstream the trachea with the development of secondary structures in the radial direction (Fig. 2).  www.nature.com/scientificreports/ The unbalance between the anterior fast and directional glottides jet and the posterior still air gives rise to double and opposites eddies. These recirculation zones were further investigated in terms of vorticity. At t = 0.2 s two opposite vortices appear near the dorsal wall of section T1 (Fig. 3), then they progressively develop until they reach a maximum value of vorticity equal to 500/s at the inspiratory peak (t = 1 s). Each vortex core is progressively shifted towards the anterior wall (T3) downstream of the trachea, while the intensity decreases. In the proximity of the 1st bifurcation, the main bronchi velocity profile (LMB, RMB) shows the peak near the inner walls, which is confirmed by other works 1,9,14,43 . Airflow pattern in downstream branches is regular and symmetric, due to the highly-laminar regime of these regions. The presence of transversal fluid velocity affect the spatial distribution of particles towards trachea wall.
Particulate phase. Steady DE (Eq. 8) of 3-25 µm particles in the oral cavity, larynx, trachea, 1st bifurcation and the four lobes are listed in Fig. 4. Increasing particle size, thus particle inertia, the filtering effect of larynx increases. DE of 7-10 µm particles is more uniform although a high concentration correctly persists in the larynx 15 (see "Supplementary Materials", Table S1). Generally, the right lung receives a greater portion of the particles, which is reasonable considering the greater flow ventilation in this region and the laminar regime simulated (< 10 L/min) 19 .
Under steady-state conditions, regional DE of 3-25 µm particles is shown in Fig. 5, respectively, for validation purposes against data in the literature.
DE of our characteristic model showed a similar trend as the empirical model of Zhou et al. 28 and other data found in the literature 3,15,[29][30][31] . Good agreement between the empirical formulation and the characteristic model data was achieved with 7 µm particles, especially at the 1st and 2nd bifurcations.
Concerning the deposition of 7 µm aerosol in a transitory simulation, total and regional DE were analysed at various injection times (Fig. 6a). In the upper regions of the model (oral cavity, larynx and trachea) DE increased with increasing injection time, reaching the maximum value at the inspiratory peak. It rapidly decreased during the deceleration phase (1-2 s) of the inhalation cycle. Particles mostly deposited at the larynx walls, where the effect of the laryngeal jet combined with geometrical features determined a peak value of DE equal to 15.07% ( t inj = 1 s). At the 2nd and 3rd bifurcations, DE reached the maximum value at t inj = 0.6 s, while injecting particles at t inj = 0.8 s seemed to be more efficient concerning the 1st (main) bifurcation. Injection at t inj = 0.8 s was the most efficient strategy also considering the amount of particles deposited in the whole model: total DE at t inj this was equal to 34.78%, which is 20.33% larger than the one observed in steady simulations.
The unsteady nature of the airflow improved TB deposition, leading to a total DE of 18.44% + 4% than the one observed under steady-conditions; detailed deposition is reported in (Fig. 6c).
The percentage of particles leaving the terminal branches at various injection times is presented in Fig. 6b. The number of particles escaped from the lower (RLL, LLL) lobes is larger than that of the particles leaving the upper and medium lobes (RUL, RML, LUL). Maximum percentages were reached at t inj = 0 − 0.2 s, except for the RML lobe, where more particles left the RML when injected at t inj = 0.8 s. On the other side, particles injected during the deceleration phase (1.4-2 s) mostly floated through the domain so that deposition or escape data were not observed.
Considering the initial homogenous distribution of particles injected at the oral inlet, Fig. 7 shows particles exiting from different lobes with different colours. Particles leaving the same lobe were closed to each other, thus forming a layered distribution which evolves in time.

Discussion
In the present work, a 3D characteristic model of the airway tree extending from the oral cavity to the 6th generation was developed to evaluate airflow dynamics and aerosol deposition under sedentary breathing physiologically conditions. Morphometric dimensions are consistent with physiological data available in the literature, extracted from in-vitro preparation of lung casts combined with CT technique. The model was validated by comparing both airflow and pressure distribution through the lobes under steady state conditions of simulation, with physiological values reported in the literature 9,16,41,42 . Steady state simulations lay the foundations of further investigations, by the implementation of a robust transitory simulation model.
Highly-unstable phenomena were also correctly captured by the model. The laryngeal jet effect propagates downstream of the trachea (Figs. 1, 2, 3) with secondary two-vortex structures, proving that the reduced-CSA glottis is the dominant site of production of separation zones 14,15,18 , then the vortex rapidly dissipates downstream the trachea down to the 1st bifurcation.
These results confirm the morphometric accuracy of our model and the importance of including anatomic details like the curvature of the upper airways, the curvatures of bronchi and their branching angles.
The unsteady model also proved to be realistic in terms of particle deposition analyses. DE of particles inhaled under sedentary breathing conditions showed good agreement with results from previous studies 15,19,[28][29][30][31] . Best fitting was provided for 7 µm particles, especially at the 1st and 2nd bifurcation (Fig. 5). Deposition in trachea shows to be less than how expected from other studies. This trend may be related to the presence of cartilaginous rings on the trachea walls, acting on the formation of boundary layer. Total deposition in transient conditions is higher than the results obtained in steady: such result points out the necessity to consider transient simulation to get a better description of the phenomenon.
Concerning the deposition of 7 µm particles, the sinusoidal airflow pattern showed a major effect on particles behaviour. The inhalation of the same particles sample at different time points led to different results in terms of both regional and total DE. The acceleration phase (0.6-1 s) proved to be the more efficient, especially on the laryngeal and 3rd bifurcations walls. Furthermore, considering the complete breathing cycle, deposition of particles inhaled in a transient regime was enhanced by 4% with respect to the constant airflow (steady) situation (Fig. 6c). Considering the injection time with the best performance, total DE was even enhanced by 20%. The improved deposition at bifurcation walls can be used to enhance the efficacy of therapeutic treatments in these regions.
Concerning particle trajectories, the distribution shown in Fig. 7 indicates that injections during the first breathing phases (0-0.2 s) drive aerosol to deep regions, so that particles leave the model. The number of particles escaped from the lower lobes (RLL, LLL) is larger than that of the particles leaving the upper and medium lobes, but there is no link between the initial position and the lobe that they escape from. Although it is impossible to predict the trajectories of particles based on their initial position, the model predicts a high percentage of 7 µm particles escaped: this aerosol is available for therapeutic deposition in deeper bronchioles of the human airways. Treatments of pathologies related to trachea are instead not optimized by this injection configuration but, based on these results, may benefit from an increase in the aerosol size.
About the modelling approach, the description of particles/droplets deposition, in conditions of laminar or quasi-transitional flow, is an issues common to a number of technical fields. Due to these heterogeneity of applications, it is not unusual to find different numerical settings, derived and dedicated to the geometrical and fluid dynamic conditions under investigation, but non-directly comparable one to each other. About this, in the field of T-junction fluid dynamics, the particles deposition on surfaces can be described by a different St formulation that represents the ratio of particle inertia to the fluid drag force 44 . Such a representation may fit better to predict particles deposition in bronchi, offering the same trend but higher values of St than the one considered in this paper. On the other hand such results would not be directly comparable to data available in the reference literature. Over the validation of a model and the comparison of new morphological and physiological conditions, further investigations on adjoining modelling solutions may consolidate the knowledge and the prediction of particles/droplets deposition in bronchi.
Our results shows the importance of including the transient behaviour of a physiologically breathing cycle to evaluate particle trajectories and deposition under sedentary conditions. Constant airflow simulations underestimate both local and total particle DE. Our model could be used to assess the efficacy of therapeutic aerosol treatment in different regions of human airways down to the 6th generation. Its applicability could be extended to different breathing conditions and to different particle types.

Conclusions
The present 3D model of the respiratory tree contains characteristic geometric details of the oral cavity, pharynx, larynx and TB tree down to the 6th generation. Morphologic and morphometric features have been carefully reproduced in the design of bronchi around the bifurcation and of the cartilaginous rings in trachea. That guarantees a physiological airflow distribution through pulmonary lobes and realistic TB pressure drops for an average adult male patient. Transient simulation of a complete breathing cycle allowed us to describe highly-unstable phenomena like the formation of a laryngeal jet and recirculation patterns and investigate the resulting aerosol deposition. Aerosol DE was analysed at various injection times and compared with value observed under constant airflow rate. The number of particles leaving each lobe was also assessed. Our results show the importance of including a transient physiologically breathing cycle on aerosol deposition analyses. Furthermore the results indicate that there not relation between the position of the inhalation across the mouth and the deposition lobe, because the vorticity in the glottides and trachea sections mix the trajectories on the transversal plane.
In conclusion, our 3D model can be used to correctly analyse airflow and aerosol behaviour for an average adult male. The proposed model was developed to evaluate airflow dynamics under sedentary breathing

Materials and methods
Anatomical model. The 3D airway model contains a characteristic geometry of the oral cavity, throat and TB tree down to the 6th generation, including physiological bifurcations and cartilaginous rings (Fig. 8a).
The reconstruction of the TB tree was based on coordinates and dimensions of the Schmidt digital model 8  where d x n P is the parent branch diameter, d x n A and d x n B are the daughter ones and x n is the exponent which determines the flow conditions in the nth bifurcation.
Section A shows the prismatic layers created near the larynx wall, Section B shows the finest tetrahedral element of the tracheal mesh. (d) Model of physiologically realistic bifurcation based on Heistracher et al. 22 ; r c is the carinal radius of curvature, R A and R B are the daughter radii and c c is the ratio of carinal curvature ratio. (CFD-Post v16, ICEM CFD v15-Ansys, www.ansys .com).
To reproduce physiologically characteristic bifurcations (Fig. 8d), every parent segment was split into two daughter tubes with a transitional zone starting from the 80% of its axial length. Then, each transition zone was smoothed defining the characteristic carinal radius of curvature r c through the equation: www.nature.com/scientificreports/ where R A and R B are the daughter radii and c c is the ratio of carinal curvature ratio equal to 0.1 23 . The tracheal section was modelled as a circular tube with a diameter equal to 1.56 cm, cut by a posterior plane at 0.65 cm from the centre 45 thus reproducing the tracheal typical C-shape due to the presence of the pars membranacea 46 . As shown in Fig. 8b, tracheal and main bronchi walls were enriched with cartilaginous rings, whose dimensions are consistent with in-vitro observations of Russo et al. 39 . The rings arrangement was made regularly perpendicular to the axial branch direction, while the rings orientation was adapted in the transitional zones to maintain their parallelism with the curvature radii. A characteristic model of extra thoracic airways was created by shaping the oral cavity, pharynx and larynx as elliptical tubes (see " Supplementary Materials", Fig. S1). Geometrical parameters of each section were extracted from Cheng's silicone model 29 . Then hydraulic diameters and cross sectional areas (CSA) were scaled in order to fit Schmidt realistic model of TB tree.
Grid generation. The 3-D characteristic airway model was divided into 24 bodies according to size and anatomical region: oral cavity and pharynx, larynx, trachea, main bronchi bifurcation (1st bif.), second (2nd bif.) and third (3rd bif.) generations. There are 64 peripheral small airways. ICEM CFD v15 software (Ansys Inc., Canonsburg, PA, USA) was employed to discretize the model with tetrahedral volume elements, triangular superficial elements and a multi-layer prism mesh to solve the high velocity gradient expected at the upper airways and larynx walls (Fig. 8c). The best size of the grid to solve this problem (trade-off between accuracy of results and calculation time) was selected by a sensitivity analysis (see " Supplementary Materials", Fig. S2) and results in a mesh of 9,826,620 tetrahedral elements and three boundary layers of hexahedral ones. Transient CFD simulations were performed to reproduce a complete sedentary breathing cycle. Such a cycle was modelled with a sinusoidal curve 4 s long, with inhalation and exhalation phases of the same duration 14 . Airflow had mean and peak values equal to 6 L/min and 9.42 L/min, respectively.
A user-defined function (UDF) was implemented to impose a time-dependent, sinusoidal velocity profile at the i-the outlet of the model, through the formula: where u i,max is the maximum velocity at the ith outflow extracted from steady-state solutions (see "Supplementary Materials"), T is the breathing period and t is time. No-slip boundary conditions and zero pressure were assumed at airway walls and at the oral inlet, respectively.
Numerical solutions and settings-particles. Particle transport equations were solved using DPM available in Fluent v16 (ANSYS, Inc., Canonsburg, PA, USA), in combination with user-defined injection patterns and boundary conditions. DPM (see "Supplementary Materials") tracks individual particles as they move through the flow using a Lagrangian approach. The Navier-Stokes equations for the continuous and discrete phases were solved separately in the steady-state condition, while they were solved together in the transient regime.
Particle deposition is usually assessed as a function of particle Stokes number St: where R 0 and U are the average radius and the mean velocity of the airflow in the parent branch, respectively and C c is the Cunningham slip correction factor defined by Hinds et al. 32 . St is typically considered to express the ratio between the particle stopping distance and a characteristic dimension of an obstacle. Thus at large Stokes numbers, particles may deviate from fluid streamlines and impact on obstacle surface, while at small numbers they tend to follow fluid streamlines 15,32 . The particles simulated in our analyses had constant diameters (3 ≤ d p ≤ 25 µm), which is the typical range of aerosol treatment. Density and viscosity were set equal to water-liquid ones, thus ρ = 998.2 kg/m 3 and viscosity µ = 0.001 kg/m s and inert behaviour was considered. C c was set equal to 1, in accordance with observations reported in Hinds et al. 32 , and Tu et al. 47 . The injection pattern consisted of 2720 particles homogeneously distributed on the oral surface, forming a disk with a diameter equal to 2.16 cm www.nature.com/scientificreports/ at a constant distance from the oral boundary. The boundary conditions for the equations governing particle motion included zero initial velocity, deposition when the particle centre reaches the wall and escape conditions at the 64 outlets. Four different diameters equal to 3, 7, 10 and 25 µm were considered to assess the effect of particle size on DE under a constant airflow rate (6 L/min). Numerical predictions of total and regional DE were compared with both empirical formulations and experimental evidences available in the literature 3,15,[29][30][31] . Tracheobronchial deposition data were validated with the empirical results of Zhou et al. 28