PHARMACEUTICAL AEROSOLS DEPOSITION DURING INHALATION , BREATH HOLDING AND EXHALATION USING CFD

In the last decades, the air pollution lead to an increase in the impact of the particle deposition in human lungs, especially with respect to respiratory diseases. Since the airway bifurcations are difficult regions to investigate, the Computational Fluid Dynamics (CFD) can be an alternative to the study of pharmaceutical aerosols that are used to the treatment of some respiratory diseases. In this study, the particle deposition was analyzed in a three-dimensional model of four ramifications or three bifurcations during three different situations: inhalation, exhalation and breath holding. The aim was to verify the medical recommendations to hold the breath during few seconds after inhale the pharmaceutical aerosols by comparison of the deposition of particles with 5 μm diameter. The numerical results of deposition during inhalation was compared to experimental data from the literature and showed good agreements. The results showed that the number of collected aerosols in the airway walls was greater for the situation of breath holding, which, in fact, confirms the medical recommendations. During the exhalation, the particles leave the domain, which would not be an interesting action after inhale the medicines.


INTRODUCTION
The particle deposition in human lungs has been largely studied due to the fast growth of the cities and the industrialization process.These events start to affect the life quality around 1970s, as stated by Brickus and Aquino Neto (1999).It is well known that the particle deposition in the human airways can be the cause of many respiratory diseases.The World Health Organization (WHO) presented the 10 leading causes of death in the world, in 2011, and three of them are indeed related to respiratory diseases (lower respiratory infections, chronic obstructive pulmonary disease and trachea, bronchus and lung cancers).These data reveal that studies on air flow and particle deposition in human lung airways need to be in constant development.Furthermore it is important to know the particle deposition pattern along the bifurcations of the human lungs since some inhalable drugs are used for the treatment of some diseases, such as bronchitis and asthma.The capability to predict the deposition of pharmaceutical aerosols on the internal surfaces of respiratory tract is important to ensure that the regions affected by the diseases receive the drugs and avoid losses of medicine aerosols.
The respiratory tract is a region with difficult access to investigate the particle deposition and aerodynamics profiles.Some in vivo experiments involve human individuals in the tests, as carried out by Bennet (1991) and Möller et al. (2009).Although these tests represent the phenomena with more precision, since it uses realistic human lung models, this investigation can be harmful to the individual's health involved in the tests.In order to avoid in vivo tests, Kim and Fisher (1999) used glass tubes to mimic some bifurcations of the human airways, by in vitro experiments.In addition to these techniques, the Computational Fluid Dynamics (CFD) can be used as in silico investigation (assisted by computational methods) since it is a tool that provides velocity, pressure and temperature fields and allows to simulate other transport phenomena, such as mass transfer between phases and chemical reactions using numerical simulation.Besides, the response of the system under different conditions, as entrance velocity and particle diameter, can be estimated without the high costs and time associated to in vitro investigation.In this context, the present study aims the analysis of particle deposition during situations of inhalation, exhalation and breath holding in a portion of human lung, using CFD.

Model Construction
A three-dimensional model of a triple bifurcation and four generations was used in this study.It is based on the Weibel (1963) convection, which identify the ramifications of the human airway by numbers, starting at the trachea, labeled by the number 0, until the last alveolus ramification, identified by number 23.In the region of the respiratory tract considered in this study, the effects of cartilaginous rings, present in the larynx, trachea and main bronchus can be disregarded.These rings can influence on air flow and trajectory of the particles, but it is difficult to measure and include these factors on the numerical simulations.The geometry of the model used in this work consists in generations 3 rd to 6 th of Weibel (1963) convention and it was constructed using ANSYS Design Modeler ® 14.0.Figure 1 shows one G3, two G4, four G5 and eight G6, completing three bifurcations.
The dimensions used for the model construction were the same used by Zhang et al. (2002) and are shown in Figure 2. A symmetry in YZ plane and x = 0 was considered and only a half of the original geometry was used to generate the mesh.Tests were performed to verify the utilization of this symmetry condition and they have been shown that the air flow fields and the particle deposition are very similar with the case in which the entire geometry was used.

Mesh Generation
The software ANSYS ICEM CFD ® was employed to discretize the domain in finite volumes, where the governing equations were solved.Hexahedral meshes with different levels of refinement were created in order to verify the grid influence over the simulation results.Hexahedral meshes are better employed in these cases since, as observed by Longest and Vinchurkar (2007), hexahedral elements are aligned to the predominant direction of the flow, avoiding the numerical diffusion commonly found when tetrahedral elements are used.For the grid independence test, it was used 5 µm diameter particles and the deposition by action of drag and gravitational forces was monitored for different number of control volumes.This test led to an ideal mesh with 305158 elements and 321564 nodes.Near the walls, the mesh was finer when compared to the rest of the computational domain to ensure the capture of all phenomena that govern the particle deposition on the internal surfaces, since in these regions the velocity and pressure gradients are higher.Some details of this finite volume mesh are shown in Figure 4.

Governing Equations
The pressure and velocity fields were obtained solving the continuity equation (Equation 1) and the momentum conservation equation (represented by Equation 2), for laminar, isothermal (310.25 K) flow.
where ρ refers to the air density .123kg m -3 ), t, to the time, v ̅, to the fluid velocity vector, p, to the pressure, ̅ to the gravity vector and ̿ , to the stress tensor, given by Equation (3).
Values of the Reynolds number of the generations showed that the flow is laminar in this portion of respiratory tract, discarding the need to use turbulence models in the simulations.The one-way coupling was used to represent the interaction between continuous (fluid) and dispersed (particles) phases, in which the solid phase has no influence on fluid phase.The diameter of the particles was 5 µm, since it is the size commonly found in pharmaceutical aerosols.The solid particles were modeled as a discrete phase using the Lagrangian approach, which allows the determination of the trajectory of the particles applying a force balance, considering drag and gravitational forces as follow: where x p , m p , d p, ρ p and v p are the position, the mass, the diameter, the density (1000 kg m -3 ) and the velocity vector of the particle, respectively, and C D is the drag coefficient given by: where a 1 , a 2 e a 3 are constants given by Morsi and Alexander (1972) and Re p is the particle Reynolds number given by:

Boundary Condition and Numerical Method
An input air flow rate of 25 L min -1 was considered at the trachea and, to obtain the flow rate and the velocity at the third generation, this flow rate was divided by a factor of 2 3 , since the airways bifurcate three times until reach G3.At inlet of the third generation it was considered that the air velocity varies with the radius according to a parabolic velocity profile.For the domain outlets (G6) a uniform pressure was assigned, and stationary no-slip condition was used for the air at the walls.Since the wall surfaces are moist and have spongy aspect, it was considered that the particles stay trapped when they collide with the walls.
A fixed time step of 2.5×10 -5 s was used for fluid and particles flow calculations.This value of time step ensures that Courant number is less than the unit, as recommended by Fortuna (2012) and Ferzi er and Perić 00 ).This dimensionless number relates the fluid velocity, the element mesh size and the time step, indicating the fluid portion that passes through the cell in the time interval considered.The convergence criterion for advancing in time was that RMS residuals were less than 10 -5 .At every time step, 53 particles were injected uniformly spaced at the surface of the third generation and the total particles released during inhalation was 2120000.It was simulated 1 second for the inhalation and, after this time, the inlet air velocity was set to zero to mimic breath holding.The trajectories of the particles that already were into the domain continue to be calculated during 2 more seconds.For the exhalation, the sixths generations are considered as entrance of the domain and the third generation is considered the outlet, following the same boundary conditions stated above.
The commercial code ANSYS Fluent ® 14.5 was employed to obtain the fluid and particles calculations.This software uses the finite volume method in the discretization of the governing equations in each control volume.The pressure-based approach was used and it creates a pressure equation derived from the continuity and momentum equations.To evaluate the diffusion and pressure gradients, the least squares cell based method was applied and, for the convective terms, the second-order difference scheme was used.Finally, the transient terms was approximated by the first order implicit formulation.The simulations were performed using 12 processing cores AMD Opteron TM 6234 and approximately 65 hours were spent to complete the calculations for the three situations (inhalation, exhalation and breath holding).

RESULTS
The results and the discussion of the simulations are shown in this section.Initially, the numerical results were compared to experimental data found in the literature.Afterwards, the air flow fields and particle deposition were evaluated for the three situations studied.

Model Validation
The experimental data used to validate the model was provided by Kim and Fisher (1999).The authors used glass tubes to resemble the generations 3 rd to 5 th , according to Weibel (1963) convention, with the same dimensions used in this work.The deposition was presented for two bifurcations separately in function of Stokes number, given by Equation ( 7): ̅ where v ̅ is the mean air velocity in the parent tube and d is the airway diameter.The Stokes number is usually employed as inertial parameter in gas-solid flows and it relates fluid and particle characteristics times.According to Crowe et al. (2012), when the Stokes number is lower than 1, the particles tend to follow fluid streamlines and, on the other hand, for Stokes greater than 1, the particles tend to deviate the fluid streamlines and, occasionally, deposit on the airway surfaces.Air flow rates of 15, 30, 60 and 96 L min -1 was used at the trachea.The particle density was 891 kg m - 3 and its diameter varied from 3 to 10 µm.These values are similar to that used by Kim and Fisher (1999).The deposition efficiencies were calculated by the ratio of the number of particles trapped in bifurcation and the particles that entered into that bifurcation.These values for the two bifurcations are shown in Figure 5 and a curve was fitted in the form of Equation ( 8 The curve for the numerical results of the each bifurcation was fitted by the least square method (r²=0,997, for the first one and r²=0,990, for the second one) and, as can be observed in Figure 5, the numerical results showed good agreement with experimental data, indicating that the model can be validated and can be used to predict the air flow and particle deposition in third to sixth generations.The deposition increase with the Stokes number, revealing that higher particle inertia leads to a greater deposition.

Air Flow Profile
The Figure 6a shows the air flow fields in plane XY calculated after 1 second inhalation.It can be observed that the parabolic profile used at the inlet of third generation leads to a greater velocity in the center of the domain.After the first bifurcation, a deflection occurs and secondary motion could be observed at fourth and fifth generations due to the pressure gradient.It can be observed that the velocity is more intense near the inner walls for the three cross-sections.For the 3-' section the secondary motion is more pronounced due to the parabolic profile used in the third generation.When a fluid flows in a curved duct, a pressure gradient appears due to the centrifugal force occasioned by the curvature.The formation of the vortex that drives the fluid from the walls to the center can be observed in Figure 6b, which shows the cross-sections velocity profiles in the three generations.This vortex appears because the fluid in the center has greater velocity than the fluid near the walls, due to the viscosity effects, and for this reason is exposed to a higher centrifugal force that drives to the walls, origins the vortical motion.

Particle Deposition
Figure 7 shows the accumulated number of deposited particles during 1 second of inhalation and 2 seconds of subsequent breath holding and exhalation.It can be observed that the deposition is greater for the breath holding than for the exhalation.This fact confirms the medical recommendations to hold the breath during few seconds after inhale the pharmaceutical aerosols.Although more particles are deposited during inspiration, it would not be interesting that the patient continue to inhale medicine aerosols, since it can lead to a waste of the drugs.The physical mechanisms that act during the inhalation period are inertial impaction and gravitational settling.When the velocity at the third generation was set to zero to mimic the breath holding, the gravitational settling was the dominant mechanism.At this moment the velocity was rapidly reduced and, consequently, the residence time became greater.The settling velocity was then rapidly reached and the particles could be deposited by action of gravitational force.

CONCLUSIONS
This study used CFD to compare the deposition of 5 µm particles in three bifurcations of human airways during inhalation, exhalation and breath periods.It was concluded that the air flow is more intense in the inner generations of model due to the parabolic profile used in the entrance of the domain.Furthermore, it was observed that the deposition is greater for the situation of breath holding than for the exhalation, which can explain the medical recommendations of hold the breath during few seconds after inhale the pharmaceutical aerosols.

Figure 1 -Figure 2 -
Figure 1 -Three-dimensional model of four generations and three bifurcations.

Figure 4 -
Figure 4 -Details of the finite volume mesh.

Figure 5 -
Figure 5 -Model validation for (a) first and (b) second bifurcation.

Figure 6 -
Figure 6 -Air flow profile (a) and cross-sections showing secondary motion (b).