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.


SM2
Figure. S1 Regions and cross sections of the respiratory system model Critical cross-sectional areas and perimeters of the upper airways model; r is the anatomical curvature radius of the throat, α is the orientation of the trachea with respect to the vertical axis (13). In addition, the trachea was oriented at an angle α=17° from the vertical axis to match CT-scan analyses by Xi and Longest and the anatomical curvature of the throat was reproduced by imposing a curvature radius r=3.27 cm 13. (CFD-Post v16 -Ansys, www.Ansys.com)

SM3
Our approach was based on the Octree-based method and Laplace smoothing to generate a high resolution surface mesh, while the volumetric elements were subsequently created with the Delaunay algorithm. Different grids consisting of 8,714,823 (Mesh 1), 9,826,620 (Mesh 2) and 10,496,004 (Mesh 3) cells were created to evaluate the grid size sensitivity for the flow solution.
A boundary layer size sensitivity was also performed to fully capture the highly complex velocity profile near the larynx-glottis walls, discretising the region at a constant distance of = 0.73 from the upper airways and larynx walls with three and ten layers of prismatic elements, respectively. Grid independency results in terms of velocity profile at the two sections A and B in the trachea. BL is the number of elements in the boundary layer. Figure S2 shows the velocity profile over the non-dimensional length in both tracheal and laryngeal sections. indicate that Mesh 1 grid was hardly sufficient to capture airflow fluid dynamics, while the increase of tetrahedral element size up to 10,496,004 cells (Mesh 3) did not change the results in terms of velocity profile in the most critical regions with respect to Mesh 2. Furthermore, three boundary layers proved to be sufficient to fully reconstruct the recirculation zones occurring at laryngeal walls. Hence, Mesh 2 with 9,826,620 tetrahedral elements and three boundary layers of hexahedral ones was adopted to perform the numerical simulations. (ICEM CFD v15 -Ansys , www.Ansys.com)

SM4 -preliminary simulations steady
Steady-state simulations were performed to evaluate airflow distribution through the model. A uniform velocity field was applied at the oral inlet with an average value of 0.2211 m/s. This value corresponds to a minute volume of 6 L/min, which is representative of sedentary breathing conditions for an average adult male (Guyton, 1977). No-slip boundary conditions were assumed at the airway walls, while constant zeropressure was prescribed at the 64 outlets of the model. Steady-state solution of the flow field was assumed convergent when the residuals of the governing equations reached values <10^(-6).

SM5
The motion of spherical liquid particles suspended in air is governed by Newton's Second Law: where , is the velocity of the i particle and the force acting on it.
Considering small particles, Reynolds numbers = | , − | µ ≪ 1, a large density ratio ≫ 1 and particle diameter d p > 1µm, most of the known particle forces other than drag force can be neglected (Zhang, 2002), thus Eq. (6) can be simplified as follows: , ⃗⃗⃗⃗⃗⃗⃗⃗ = (⃗⃗⃗ − , ⃗⃗⃗⃗⃗⃗ ) (7) where , and are the components of particle and fluid velocity in the − ℎ position, respectively. For spherical particles, the drag force per unit mass is defined as: where C D is the drag coefficient defined by Morsi and Alexander (Morsi, 1972), d p is the particle diameter and C c is the Cunningham slip correction factor defined by Hinds et al. (Hinds, 1999). With particle Reynolds numbers Re p < 1 and very small Stokes numbers (St < 0.25), several authors found that particle motion can be basically achieved with a first-order correction to fluid element motion equation (Comer 2000) (Zhang 2001Zhang 2002). For example, Comer observed that the largest influence of the Cunningham slip correction factor was a 2% change in particle DE at St = 0.24 (Comer, PhD thesis publication). In this sense, is used for particles in the nanometric range, where the primary phase cannot be considered continuous (Tu, 2013). Furthermore, the correction to Stokes law for non-rigid spheres as water droplets is generally insignificant (Hinds 1999; Tu 2013; Fluent Guide, Fluent 16, Ansys Inc.).  Table S1. Deposition efficiency (DE) of 3-25µm particles injected under steady conditions (Q=6 L/min). RUL and RLL are the right upper and lower lobes, LUL and LLL are the left upper and lower ones. Grey cells represent the maximum DE of each region of the model, bold font highlights the maximum DE for various particle diameters.