of Chemical Engineering A CFD STUDY OF DEPOSITION OF PHARMACEUTICAL AEROSOLS UNDER DIFFERENT RESPIRATORY CONDITIONS

Respiratory diseases have received increasing attention in recent decades. Airway bifurcations are difficult regions to study, and computational fluid dynamics (CFD) offers an alternative way of evaluating the behavior of pharmaceutical aerosols used in the treatment of respiratory disorders. In this work, particle deposition was analyzed using a three-dimensional model with four ramifications (three bifurcations), under different respiratory conditions: inhalation, exhalation, and breath holding. The main aim of the work was to verify the medical recommendation to hold one’s breath during a few seconds after inhaling pharmaceutical aerosols, rather than exhaling immediately after the inhalation. The deposition of particles with 5 μm diameter was considered. The results showed that the number of aerosols collected on the airway walls was higher for the situation of breath holding, which supported the medical recommendation.


INTRODUCTION
Studies of particle deposition in the human lung have largely been developed as a response to air pollution and associated effects on the quality of life, resulting from the rapid growth of urban areas and industrialization (Brickus and Aquino Neto, 1999).The deposition of particles in human airways can be the cause of many respiratory diseases.According to the World Health Organization (WHO), in 2013 three out of the ten leading causes of death worldwide were related to respiratory diseases (lower respiratory system infections, chronic obstructive pulmonary disease, and tracheal, bronchial, and lung cancers).These data support the need for ongoing studies of air flow and particle deposition in human lung airways.Furthermore, knowledge of the pattern of particle deposition along the bifurcations of the human lung is indispensable, because some inhalable drugs are used for the treatment of respiratory diseases such as bronchitis and asthma.The capability to predict the deposition pattern of pharmaceutical aerosols on the internal surfaces of the respiratory tract is important in order to ensure that the regions affected by the diseases receive the drugs, as well as to avoid losses of medicinal aerosols.
The respiratory tract is a region that is difficult to access for investigation of particle aerodynamics and deposition patterns.Examples of in vivo studies include the work of Bennet (1991) and Möller et al. (2009).Although such investigations can describe the phenomena with high precision (since real human lungs are used), the techniques employed can be potentially harmful to the health of the individuals involved in the tests.In order to avoid in vivo tests, Kim and Fisher (1999)

used in vitro experiments
Brazilian Journal of Chemical Engineering with glass tubes to mimic some of the bifurcations of the human airways.In addition to these methods, computational fluid dynamics (CFD) can be used in in silico investigations employing computational techniques.CFD is a tool that enables the simulation of velocity, pressure, and temperature fields, as well as other transport phenomena.Comparison of the numerical results with experimental data can be used to evaluate the reliability of the mathematical model and its ability to provide a realistic representation of the physical phenomena.In addition, the response of the modeled system under different conditions, considering parameters such as entrance velocity and particle diameter, can be determined without the high costs and lengthy procedures associated with in vitro investigations.
Numerical studies of the fluid dynamics and particle deposition inside the lungs first requires the construction of a morphological model representing the respiratory tract.The scheme that is most widely used to describe the human airways was created by Weibel (1963).This identifies the lung ramifications (called generations) using numbers, starting at the trachea, labeled with the number 0, and continuing until the final alveolar ramifications, identified by the number 23.These geometries can be generated using real images, enabling the creation of more realistic models (Takano et al., 2006;Xi and Longest, 2007;Luo and Liu, 2008;Imai et al., 2012).Although this technique can reproduce the real dimensions of the human airways (taking into consideration the precision of the camera), the imperfections of the ramifications can complicate the numerical mesh generation and, consequently, can cause numerical errors to occur.This is the main reason to treat the bifurcated ducts of the geometric model as cylinders.The average dimensions of each generation (for middle-aged men and women) can be found in ICRP (1994).
The literature contains a substantial body of information concerning pharmaceutical aerosol deposition.Longest and Holbrook (2012) presented a brief history of the use of computational methods to investigate drug delivery in human airways.Tena and Clarà (2012) presented a summary of the main parameters affecting the deposition of particles in human lungs, and provided details of drug deposition for the treatment of respiratory diseases, including a discussion of devices used for the administration of inhaled drugs.Yang et al. (2014) published a review that focused on inhaled dry powder aerosols and discussed the respiratory and particle parameters that most affected the operation of the devices used.Some studies have found that computational methods can help in the optimization of nebulizers and inhalers (Longest and Hindle, 2009;Xie et al., 2010;Cui et al., 2014).Chen et al. (2012a) studied the gassolid flow in an obstructed airway, a condition commonly found in patients with COPD (chronic obstructive pulmonary disease).It was found that there was higher recirculation and secondary motion of the air and particle flow in the obstructed airways and in the region below this generation, compared to the flow into healthy lungs.This condition can significantly increase drug delivery to the wall surfaces.Similar studies have involved stochastic models of ramifications affected by emphysema (Sturm and Hofmann, 2004), where it was found that deposition was reduced in the affected airways.
Most of the studies undertaken to understand the influence of the particle properties, respiratory action, and morphological characteristics have considered the inhalation step (Liu et al., 2003;Kleinstreuer et al., 2007;Longest and Vinchurkar, 2007b;Chen et al., 2012a;Chen et al., 2012b;Saber and Heydari, 2012;Piglione et al., 2012) or the exhalation stage (Balásházy et al., 2002;Zhang et al., 2002;Nowak et al., 2003;Longest and Vinchurkar, 2009).Few investigations have considered a breath holding step as a strategy to increase the deposition of inhaled drug aerosols.Inthavong et al. (2010) conducted a numerical study of the deposition of pharmaceutical aerosols in the first six generations of the lung, considering breath holding.However, the velocity fields and the influence of this action on particle deposition were not explored in detail.Imai et al. (2012) also studied the deposition of particles during breath holding and observed an increase of the deposition of 5 µm particles, compared to the inhalation step.The effect of breath holding was not compared with the exhalation step in the studies of Inthavong et al. (2010) and Imai et al. (2012).
In the context of pharmaceutical aerosols delivered to the respiratory tract, the contribution of the present study is a CFD analysis of particle deposition in a portion of the human lung during situations of inhalation, exhalation, and breath holding.Drug deposition during these three situations was compared in order to verify the medical recommendation to hold one's breath after the inhalation of aerosols from nebulizers and inhalers.

MATHEMATICAL MODEL
The continuity and the momentum conservation equations, represented by Equation (1) and Equation (2), respectively, were solved for a three-dimensional gas flow, employing an Eulerian approach.The Reynolds number for all the four generations was calculated based on the mean velocity of the gas and the duct diameter in each bifurcation and for each case studied.The results showed values ranging from 173 to 669, confirming laminar flow in this portion of the respiratory tract.Thus, turbulence models were not required in the simulations.
In Equations ( 1) and ( 2), ρf is the air density, t is the time, v  is the velocity vector, p is the pressure, g  is the gravity vector, and    is the stress tensor, given by Equation (3): where µ is the air viscosity and I   is the unit tensor.
The values of the physical properties of the fluid were chosen based on the normal temperature of the human body (approximately 310.15 K).The values of these parameters are presented in Table 1.The drug particles were modeled as a discrete phase using a Lagrangian approach, which allows determination of the trajectories of the individual particles by employing a force balance.Since the particulates occupy less than 1% of the volumetric fraction, when compared to the continuous phase, it is reasonable to assume that one-way coupling can represent the interaction between the continuous (fluid) and dispersed (particles) phases.In this approach, only the fluid has an influence on the trajectories of the particles, and the discrete phase has no influence on the fluid calculations.The ratio between the momentum response time of a particle (τV) and the time between collisions (τC) gives an indication of whether the flow is dense or dilute (Crowe et al., 2012).If this ratio is less than one, the flow is considered to be dilute and fluid forces control the particle motion.These characteristic times can be calculated as: where n is the number density of particles, dp is the particle diameter, vr is the relative velocity of a particle with respect to the others, and ρp is the particle density.Considering the values for the physical properties of the particle and fluid phases shown in Table 1, the value obtained for the ratio was 0.0031, which confirmed that the use of one-way coupling was suitable in the present work.The particle diameter used was 5 µm, because this is the size commonly found for pharmaceutical aerosols (Amaral, 2010).The physical properties of the discrete phase are shown in Table 1.The density considered for the particles was approximately the same as found for drug solutions.The force balance used to estimate the trajectory of the particles included drag and gravitational forces, as follows: where p v  and Rep are the velocity vector and the dimensionless Reynolds number of the particle, respectively, and CD is the drag coefficient, given by: where a1, a2, and a3 are constants described by Morsi and Alexander (1972), and Rep is the particle Reynolds number, defined by: According to Finlay (2001), the ratio between the distance that a particle diffuses (xd) and the distance that a particle settles (xs) is a measure of the importance of Brownian diffusion relative to sedimentation.When this ratio is less than 0.1, diffusion can be neglected, for a particle of size dp.This ratio can be calculated by the equation: where kB is the Boltzmann constant, T is the absolute temperature, and CC is the Cunningham correction factor.For the particular case of this paper, the value of the ratio was 0.025, indicating that Brownian motion could be neglected in calculation of the particle trajectories.

SIMULATIONS
The equations described in the previous section for the mathematical modeling were numerically solved using the commercial code ANSYS Fluent 14.5 to obtain the fluid and particle profiles.To solve the partial differential equations, this software applies the finite volume method, where the computational domain is discretized into finite control volumes for which the governing equations are solved and the variables are calculated.In the case of this study, the continuity and momentum equations were manipulated to create a derived pressure equation, using the pressure-based approach.To evaluate the diffusion and pressure gradients, the least squares cell-based method was employed to interpolate the diffusion and pressure gradients, while the convective terms were interpolated by the second order upwind difference scheme.Finally, the transient terms were approximated by the first order implicit formulation.The system of algebraic equations was solved by a segregated algorithm available in ANSYS Fluent named SIMPLE (Semi Implicit Method for Pressure Linked Equations).Double precision was used in all calculations for the set for simulations.

Boundary Conditions and Simulation Details
The trachea input air flow rate was 25 L min -1 , as described by ICRP (1994).The airways bifurcated three times until reaching G3, so the flow rate at the third generation could be obtained by dividing the flow rate at the trachea by a factor of 2 3 .The air flow rate and the diameter of the ramification were then used to calculate the velocity at the third generation.The variation of the air velocity with the radius was also considered, using a parabolic velocity profile where the velocity was highest at the center of the respiratory duct and lowest at the walls.This was achieved with a user-defined function (UDF) implemented in the software to consider the variation of the inlet velocity at G3.A zero relative pressure was adopted for the domain outlets (G6).Finally, a stationary no-slip condition was used for the air at the walls.Moreover, since lung wall surfaces are moist and have a spongy texture, it was assumed that the particles remained trapped after collision with the walls.
A fixed time step of 2.5×10 -5 s was used for the fluid and particle flow calculations, which ensured that the Courant number was less than unity, as recommended by Fortuna (2012) and Ferziger and Perić (2002).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 the RMS residuals were less than 10 -5 .
At every time step, 53 particles were injected (uniformly spaced) at the inlet surface of the third generation, and the total number of particles released during inhalation was 2,120,000.The particle trajectories were calculated by the Discrete Phase Model (DPM) available in the software.In this method, the particles are considered as points and are tracked individually.The simulation used one second for the inhalation, after which the inlet air velocity was set to zero to mimic breath holding.Calculation of the trajectories of the particles that were already in the domain was continued during two more seconds.For the exhalation, the sixth generations were considered as the entrance of the domain, and the third generation was considered the outlet, following the same boundary conditions described above.Similarly, two seconds were used in the simulation of the exhalation step in order to obtain the trajectories of the particles that were already inside the domain at the end of the inhalation.The simulations were performed using 12 processing cores of an AMD Opteron TM 6234, and approximately 65 hours were required to complete the calculations for the three situations (inhalation, exhalation, and breath holding).

Model Validation
The experimental data used to validate the model were obtained from Kim and Fisher (1999) The curves for the numerical results obtained for each bifurcation were fitted by the least squares method, and r² values of 0.997 and 0.990 were obtained for the first and second bifurcations, respectively.The numerical results showed good agreement with the experimental data (Figure 5), validating the model for use in prediction of the air flow and particle deposition in the third to sixth generations.Furthermore, the deposition increased with the Stokes number, confirming that higher particle inertia resulted in greater deposition.

RESULTS
The air flow field calculated after 1 second of inhalation is shown in Figure 6   Figure 9 shows the accumulated number of particles deposited during 1 second of inhalation and 2 seconds of subsequent breath holding or exhalation, from which it can be seen that there was greater deposition for breath holding than for exhalation.This finding supports medical recommendations to hold the breath for a few seconds after inhaling pharmaceutical aerosols.If the patient continued the breathing cycle of inhalation followed by exhalation, drug accumulation would be less efficient because the particles would pass through the airway bifurcations.Although more particles are deposited during inspiration, compared to breath holding and exhalation, there is a significant increase in deposition when breath holding is performed.

A
three-dimensional model of a triple bifurcation and four generations was constructed in order to represent a portion of the bronchial region of the respiratory tract, based on the Weibel (1963) convention.In this study, the bifurcation ducts corresponded to the third to sixth generations.Kim and Fisher (1999) provided data obtained in experimental tests using glass tubes with the same dimensions as the generations considered in this study, which enabled validation of the numerical results obtained in the simulations performed in the present work.The geometry of the model was constructed using ANSYS Design Modeler ® 14.5 and the final model is shown in Figure 1.The dimensions used for the model construction were the same as those employed by Zhang et al. (2002) and are shown in Figure 2.

Figure 1 :
Figure 1: Three-dimensional model with four generations and three bifurcations.

Figure 2 :
Figure 2: Dimensions of the model (in centimeters).The physical phenomena and the geometry are symmetric with respect to the YZ plane, so only half of the original model was used to generate the mesh.Simulations were performed to check this condition of symmetry, considering the full and symmetric models under the same boundary conditions.The deposition efficiencies showed a difference of 2%, demonstrating that the symmetry condition could be used without any loss of quality of the results.This strategy is useful because fewer elements are required to compose the mesh, hence decreasing the time required to complete the simulations.The ANSYS ICEM CFD ® 14.5 software package was used to discretize the domain into finite volumes for which the governing equations were solved.Hexahedral meshes with different levels of refinement were created in order to determine the influence of the grid on the simulation results.According to Longest and Vinchurkar (2007a), hexahedral meshes

Figure 3 :
Figure 3: Results of grid independence tests.

Figure 4 :
Figure 4: Details of the finite volume mesh.
v is the mean air velocity in the parent tube, which is the main duct that gives rise to the other two in a bifurcation, and d is the airway diameter.The mean velocity was calculated considering the air flow rate and the cross-sectional area at the parent duct.The Stokes number is usually employed as an inertial parameter that describes the fluid and particle characteristic times in gas-solid flows.According to Crowe et al. (2012), when the Stokes number is lower than 1, the particles tend to follow the fluid streamlines, while for Stokes numbers greater than 1, the particles tend to deviate from the fluid streamlines and can deposit onto the airway surfaces.The model validation simulations employed the parameters and boundary conditions described by Kim and Fisher (1999), as shown in Table 2.The deposition efficiencies were calculated as the ratio of the number of particles trapped in a bifurcation to the number of particles that entered the bifurcation.The numerical and experimental values obtained for the two bifurcations are shown in Figure 5, and a curve was fitted in the form of Equation (11), where a and b are constants of the model (Kim and Fisher, 1999):

Figure 5 :
Figure 5: Model validation for the first (a) and second (b) bifurcations.
. The parabolic profile used as the boundary condition at the inlet of the third generation resulted in a higher velocity in the center of the computational domain following the entrance.After the first bifurcation, a deflection occurred and secondary motion could be observed at the fourth and fifth generations, due to the pressure gradient.When a fluid flows in a curved duct, a pressure gradient arises due to the centrifugal force caused by the curvature.The formation of a vortex driving the fluid from the walls to the center can be observed in the cross-sectional velocity profiles in the three generations (Figures 6(b) to 6(d)).This vortex appears because the fluid in the center has greater velocity than the fluid near the walls, due to the viscosity effects.The central fluid portion is therefore exposed to a higher centrifugal force that drives towards the walls, giving rise to the vortex motion.Consequently, the velocity was more intense near the inner walls, for all three cross-sections.The secondary motion was more pronounced for the 3-3' section, due to the parabolic profile used in the third generation.The velocity profile at the beginning of the breath holding is shown in Figure 7.It can be observed that the velocity decreased rapidly across the entire Brazilian Journal of Chemical Engineering Vol. 33, No. 03, pp.549 -558, July -September, 2016 domain when the inlet velocity was set to zero in order to mimic breath holding.

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

Figure 7 :
Figure 7: Air flow profile during the beginning of breath holding.The accumulated number of particles deposited over time during the inhalation situation is shown in Figure 8.A linear profile was obtained, which was the main reason for the choice of the simulation time.Medical recommendations often suggest that the patient should inhale a drug during 5 to 10 seconds.Here, the linear profile observed for simulation during 1 second indicates that the results could be extrapolated for the next seconds, allowing particle deposition to be determined for any duration of inhalation.Figure9shows the accumulated number of particles deposited during 1 second of inhalation and 2 seconds of subsequent breath holding or exhalation, from which it can be seen that there was greater deposition for breath holding than for exhalation.This finding supports medical recommendations to hold the breath for a few seconds after inhaling pharmaceutical aerosols.If the patient continued the

Figure 8 :
Figure 8: Number of particles deposited during the inhalation stage.

Figure 9 :
Figure 9: Number of particles deposited during inhalation, breath holding, and exhalation.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, in order to mimic breath holding, the dominant mechanism became gravitational settling.At this moment, the velocity decreased rapidly and the residence time consequently became longer.The settling velocity was then rapidly reached and the particles could be deposited by the action of gravitational force.The Stokes number, described above, can also be defined by the ratio between the time that a particle takes to respond to the modification in the fluid velocity (τv) to a characteristic time of the fluid (Crowe et al., 2012): = v v St d 

Figure
Figure 10: Streamlines after 1 second of inhalation (a) and after 0.5 of exhalation (b) and breath holding (c).Comparison of the fluid streamlines during exhalation (Figure 10(b)) and breath holding (Figure 10(c)) indicates that, in the first case, the particles probably flow out of the domain without any significant contact with the walls, because they tend to follow the fluid streamlines.On the other hand, during breath holding, some vortices can be observed and the particles are likely to recirculate in the computational