Importance of Physical and Physiological Parameters in Simulated Particle Transport in the Alveolar Zone of the Human Lung

The trajectory and deposition efficiency of micron-sized (1–5 μm) particles, inhaled into the pulmonary system, are accurately determined with the aid of a newly developed model and modified simulation techniques. This alveolar model, which has a simple but physiologically appropriate geometry, and the utilized fluid structure interaction (FSI) methods permit the precise simulation of tissue wall deformation and particle fluid interactions. The relation between tissue movement and airflow in the alveolated duct is solved by a two-way fluid structure interaction simulation technique, using ANSYS Workbench (Release 16.0, ANSYS INC., Pittsburgh, PA, USA, 2015). The dynamic transport of particles and their deposition are investigated as a function of aerodynamic particle size, tissue visco-elasticity, tidal breathing period, gravity orientation and particle–fluid interactions. It is found that the fluid flows and streamlines differ between the present flexible model and rigid models, and the two-way coupling particle trajectories vary relative to one-way particle coupling. In addition, the results indicate that modelling the two-way coupling particle system is important because the two-way discrete phase method (DPM) approach despite its complexity provides more extensive particle interactions and is more reliable than transport results from the one-way DPM approach. The substantial difference between the results of the two approaches is likely due to particle–fluid interactions, which re-suspend the sediment particles in the airway stream and hence pass from the current generation.


Introduction
Nowadays, for chronic obstructive and especially asthmatic diseases in the pulmonary respiratory system, anti-inflammatory medications are used in the form of inhalation drugs.Because of their high effectiveness and low side effects, inhaler drugs taken via respiratory airways are the main approach for the treatment of lung diseases.
Many mathematical models and computation techniques have been developed or modified for predicting accurately the transport of particles in the respiratory tract.Various modelling approaches have been applied for analysing the aerosol deposition performance within the pulmonary system [1,2], for example semi-empirical [3], symmetric generation [4], trumpet [5], asymmetric multiple path modes [6], stochastic, asymmetric generation [2] and computational fluid dynamics (CFD) models [7].Unfortunately, the results of these models vary primarily due to different lung inlet morphometric conditions such as breathing periods, pressure gradients and mathematical modelling techniques.Some researchers have investigated the particle deposition efficiency in deformable domains by determining the pre-defined wall motions of the lung.For example, Henry et al. applied sinusoidal oscillations for a fully-alveolated alveolar duct [8], but they neglected Brownian diffusion forces.However, some studies have accounted for the Brownian force and compared the root mean square (RMS) displacement of a Brownian particle to the theoretical diffusion predicted by Einstein theory for validating a particle motion algorithm [9].The gravitational deposition of particles in rhythmically expanding and contracting alveolar models has been analysed and it was shown that wall motion is a cardinal factor in particle transport [10].In this regard, fluid structure interaction models can be developed for more accurate results and to investigate how breathing pattern models and lung tissue mechanical properties affect deep-lung flow fields and particle dynamical transport [9].
In this article, we examine computationally particle transport in an idealized human alveolar model.The model has several features and advantages such as accurate simulation of tissue wall deformation and particle-fluid interaction collisions.Although the developed alveolar model used in this study has a simple geometry, it is demonstrated to be physiologically-relevant by utilizing the actual relationship between tissue motion and airflow in the alveolar.This relationship is solved using the two-way fluid structure interaction (FSI) simulation technique.The main parameters affecting particle transport such as particle diameter, tissue visco-elasticity and tidal breathing period are investigated with the aid of CFD analyses for the developed model.In addition, the fluid flows and streamlines are calculated for the present flexible model and the two-way coupling particle trajectories are investigated relative to one-way particle coupling systems.The results are applicable to pharmacology applications and research.

Model Development and FSI Approach
The airspace dimensions have been derived from Haefeli-Bleuer and Weibel morphometric data for human pulmonary acinus and are representative of a total lung volume of 3.5 L. The actual dimensions of first for generations of the asinus model (G18-22) are calculated to be 265 µm for the lumen diameter, 575 µm for the outer diameter, 600 µm for the duct length and 150 µm for the alveolar opening [11].Haefeli-Bleuer and Weibel proposed an idealized model of the alveolated lung in which they gave the total surface area of alveolar walls in each generation per acinus.We set the number of axial and radial alveolar walls to get the closest match between the total alveolar surface area of their model while keeping the alveolar width similar to the alveolar depth.Our model assumes that all generations have identical linear dimensions due to the fact that the effort required in meshing the bifurcation is large.This assumption is a valid one due to the similarity in linear dimensions over 18-22 generations measured by Haefeli-Bleuer and Weibel.
In previous studies on the alveolar zone of the human lung, the alveoli walls have been assumed rigid.However, there are significant differences between real (in vivo) and computed (in silico) results because alveoli have visco-elastic tissue.The model developed here, shown in Figure 1, is three-dimensional, fully alveolated and has a realistic alveolar geometry consisting of radial alveolar walls.This type of wall has generally not been considered in previous models [12].
In addition, tissue tensions are accounted for here in modelling the lung.The lung parenchyma transfers the tissue tensions that produce a pressure gradient in the surrounding alveoli.This load expands the alveoli and establishes a sub-ambient pressure (vacuum) within the lungs.Only a few studies have employed visco-elastic walls and calculated tissue stresses, but these studies were performed only for two-dimensional models [9].Moreover, as it is important to choose suitable physiological parameter values such as magnitude of the pressure in the modelling, we utilize the following equation for sub-ambient pressure (vacuum), which was reported by Dailey and Ghadiali [9]: where ∆P = 2000 dyn/cm 2 is pressure gradient, and t is time and λ TB is the breathing period as well as the baseline values for the tissue properties from the mentioned study [9].
Appl.Sci.2017, 7, 113 3 of 22 following equation for sub-ambient pressure (vacuum), which was reported by Dailey and Ghadiali [9]: where ∆P = 2000 dyn/cm 2 is pressure gradient, and t is time and λ TB is the breathing period as well as the baseline values for the tissue properties from the mentioned study [9].

Fluid-Structure Interaction Equations
The ANSYS finite element (FE) package is used to solve the fluid structure interaction simulation.While the simulation is in process, the solver simultaneously uses a mixed discretization with Lagrangian-Eulerian (ALE) formulas.The Lagrangian and Eulerian formulations have been used for deformable walls and fluid domains, respectively.Details of the FSI formulations are available elsewhere [11].The visco-elastic model draws on the Kelvin-Voigt visco-elastic model shown in Figure 2.

Fluid-Structure Interaction Equations
The ANSYS finite element (FE) package is used to solve the fluid structure interaction simulation.While the simulation is in process, the solver simultaneously uses a mixed discretization with Lagrangian-Eulerian (ALE) formulas.The Lagrangian and Eulerian formulations have been used for deformable walls and fluid domains, respectively.Details of the FSI formulations are available elsewhere [11].The visco-elastic model draws on the Kelvin-Voigt visco-elastic model shown in Figure 2. following equation for sub-ambient pressure (vacuum), which was reported by Dailey and Ghadiali [9]: where ∆P = 2000 dyn/cm 2 is pressure gradient, and t is time and λ TB is the breathing period as well as the baseline values for the tissue properties from the mentioned study [9].

Fluid-Structure Interaction Equations
The ANSYS finite element (FE) package is used to solve the fluid structure interaction simulation.While the simulation is in process, the solver simultaneously uses a mixed discretization with Lagrangian-Eulerian (ALE) formulas.The Lagrangian and Eulerian formulations have been used for deformable walls and fluid domains, respectively.Details of the FSI formulations are available elsewhere [11].The visco-elastic model draws on the Kelvin-Voigt visco-elastic model shown in Figure 2.  The Kelvin-Voigt model has both elasticity and viscosity properties simultaneously.A purely viscous damper and a purely elastic spring are connected in parallel in the Kelvin-Voigt model (see Figure 2).
σ Total = σ D + σ S (2) In Equation ( 4), η is the viscosity.The boundary conditions at the interface regions between the two domains must be satisfied by coupling the fluid and structure domains at FSI modeling technique as shown below.
In Equation ( 4) the fluid and structure nodal displacements are shown by d f i and d s i , respectively.In addition, n j is the interface normal vector and σ ij is the Kronecker delta function.The fluid and solid stress tensors are shown by σ s ij and σ ij , respectively.The fluid velocity continuity condition in Equation ( 9) needs that the structure and fluid domains have equal nodal velocities at the movable fluid-structure boundary.The boundary conditions at the interface regions of fluid and structure domains must be satisfied by coupling the fluid and structure domains with the FSI technique [11].For this, we set up a pair of coupled systems consisting of a Transient Structural system and a Fluid Flow (CFX) system for performing the two-way FSI analysis in ANSYS Workbench 16.In this research, we utilize an independent discretization model of the fluid and structure domains.Our structure meshes are suppressed during the CFX system set up by Transient Structural system, and consists of multi-block and body-fitted hexahedral elements.While the fluid meshes set up by the CFX system contain 8304 total elements and 4561 nodes, there are 7064 total elements and 25,341 nodes in the solid model.The FSI technique defines the interface between the structure domain in the Transient Structural system and the fluid domain in the CFX system.Data of solid and fluid conditions are interchanged across this interface during the implementation of the simulation.
To ensure valid results, mesh independency for the obtained results was investigated by resolving the airflow field for a finer mesh with 1,104,000 cells for the 18th generation.There was no noticeable difference between the fine and course meshes in terms of flow field and deposition of 1-5 µm diameter particles.Specifically, the results were observed to vary approximately 0.7% and 0.2% for 1 and 5 µm diameter particles, respectively.

Flow and Particle Transport Simulation
The airflow field was developed by a CFD code using the ANSYS Workbench finite element algorithm.In this study, we used the laminar, incompressible and isothermal Navier-Stokes equations with air specifications at 37 • C for modelling the airflow conditions in the alveolar region.Because of the low Reynolds number and laminar flow regime, the velocity profile at the inlet of the model was assumed to be parabolic, and a no-slip boundary condition was applied at the wall boundaries.
The volumetric flow rate .Q G for various generations was scaled from the tracheal flow rate .Q T = 500 mL/s as follows: .
where G is the generation number.The application of Lagrangian tracking in the modelling of particle transport involves the integration of individual particles that are tracked from their injection point until they escape the domain.In addition, just as the fluid affects particle transport through forces, there is a counteracting influence of the particle on the fluid, i.e., there is a coupling between the fluid and the particles.One-way coupling occurs if the fluid is permitted to affect particle dynamics but the particles do not influence the fluid pathways, while two-way coupling exists if the particles influence the fluid behaviour reciprocally.
The particle-tracking model implemented in this study, which is referred to as the one-and two-way particle-fluid discrete element method (DEM), uses the Euler-Lagrange approach.Therefore, the mutual effect of airflow and particles is calculated along with the mutual interactions between spherical particles by means of binary collisions [13].The Lagrangian approach is used for tracking the particles that are computed with the post-processor of the CFD package.In order to compute the particle trajectories, we use the Langevin equation, which balances the mass acceleration of the particle with the forces acting on it if the particulate density is large in comparison to the fluid density: where u p , F D , F Br , mg and F VM denote velocity of the particles, drag force, force resulting from Brownian motion, gravitational force and virtual (or added) force, respectively.Moreover, the inter-particle contact forces resulting from the DEM approach are determined by the F C,ij force.The governing equations for forces are described elsewhere [14][15][16][17].
where u, and A p are fluid velocity, fluid density and cross section area of the particles and for Re p ≤ 1000: The C D is the drag coefficient.Re p is particle Reynolds number and is expressed as follows: In this equation, d p is the particle diameter and µ is air viscosity.The Brownian term F B i defined by Equation ( 15) is stochastic random vector and ∆t is the solution process time step.
where S 0 , ζ i and v are the amplitude of the white noise process, zero-mean unit-variance Gaussian random vector and the kinematic viscosity of the air, respectively.In addition, k = 1.38 × 10 −16 is the Boltzmann constant and T = 37 • C is the temperature of the air.The virtual mass force is shown by F VM in Langevin equation.This force is very effective when particles fill the major part of the volume of the fluid.
In calculating this force, the C VM and M F coefficients are normally set to 0.5 and 1, respectively.The inter-particle contact forces are divided to normal contact and tangential contact forces.
The normal contact force was computed by using the Hertz-Mindlin no-slip model [5].The contact force can be identified by describing particle contacts as damped harmonic oscillators.
where δ n , E eq , R eq are normal displacement, equivalent Young's modulus and equivalent radius, respectively.
where E, ν and R are Young's modulus, Poisson ratio and radius for elements, respectively.In addition, U n rel is the normal component of the relative velocity, e is the coefficient of restitution, m eq is the equivalent mass and S n is the normal stiffness defined as: S n = 2E eq δ n R eq (23) The m is mass of particle elements.The tangential force is defined by tangential elastic force and tangential damping force.
where U t rel is the tangential component of the relative velocity, δ t is the tangential displacement, µ p is the coefficient of static friction and S t is the tangential stiffness which is given as: In this equation the G eq is the equivalent shear modulus which is defined as

Boundary Condition
To avoid entrance effects for numerical model validation, a parabolic velocity was employed at the inlet, which was created via program.The averaged inlet velocity was determined from inhalation flow rates at the mouth.
The particle-wall interaction boundary condition was assumed to be a "100% trapped wall", because of the existence of mucus layers which coat the inner wall of lung airways.In addition, uniform (zero) gauge pressure was applied at the terminal outlets and also a no-slip boundary condition was applied at the wall boundaries.
To conserve mass flow during simulation of breathing, we need to introduce same amount of mass flow left integrated over time from the outlet.In order to accomplish this, we require tracking amount of mass flow leaving alveolated duct and integrate it over the tracking time.In our calculations we have automatized this process by computing mass flow leaving the domain from each of the outlet and we automatically switch boundary conditions and apply the same flow rate to conserve mass flow.We were able to successfully simulate this for four-generations (G18-22).

Analysis Description
The behaviour of micron-sized particles in the human alveolar is investigated here considering particle-fluid collisions, accounting for interactions between small particles and the airflow field.Initially the full FSI simulation was developed for modelling the transient flow fields for four breathing periods.A mixed discretization with Lagrangian-Eulerian formulas was utilized for setting up the FSI simulations in the ANSYS FE package.We also demonstrated the importance of effective parameters (such as particle diameter, tissue visco-elasticity and tidal breathing period) in simulated particle transport in the alveolar zone of the human lung.For the simulations, we used seven tissue models (using the Kelvin Voigt visco-elastic model with four appropriate values for a linear elastic tissue in the range of 20,000 < E < 50,000 dyn/cm 2 ; and three appropriate values for visco-elastic tissue 0 < η < 60,000 g/cms) and four breathing periods (for standard and common treatment groups λ TB = 5, 7.5, 11.25, 15 s).
The CFX system was used for creating the fluid domains including one-way coupling and two-way (full) coupling between the continuous phase and particles.For calculating the influence of particle phase on the continuous airflow pathways we ultimately need to use the full coupling system.However, the full coupling system has higher CPU cost than one-way coupling.Therefore, to solve these problems, two sets of identical particles including one-and two-way coupling were created separately for the optimization of CPU usage.The particle number rate parameter for applications involving tracking of discrete particles, it is not practical to track all physically existing particles.Instead representative particles, or parcels, are used to track these discrete particles.Each representative particle characterizes a certain number of actual particles.The actual number of particles represented by the representative particle is called the particle number rate.It is determined from the mass flow rate assigned to the representative particle divided by the mass of an actual particle.In addition, the minimum number rate that provides stable number rate (N) is defined for two-way coupling particles.
The 1-5 µm diameter particles, which have a density of 1 g/cm 3 , were tracked from inlet to outlet of the alveolated duct model during the mentioned inspiration times.Finally, we assessed the time-dependent flow data and particle-tracking results achieved via the ANSYS Workbench program.The detailed mathematical model used in this study can be obtained from the corresponding author at the above email address.

Results and Discussion
The particles injected in the developed model were distributed uniformly at the inlet region.As shown in Figure 3, the critical injected particle number rate N was stabilized by the particle deposition percentage, which must be independent of the particle number rate.This modification was carried out for one-and two-way coupling particles.The numerical simulations were performed on a local Dell Precision system with 12 GB RAM and four 3.33 GHz CPUs.
When the fluctuations of deposition became stable for a minimum value of N, we adjusted the particle number rate to this value.When N reaches 1500 and 4000 for one-and two-way coupling, respectively, the results are shown in Figure 3 to be approximately independent of N.

Flow Field and Particle Trajectories
For the rigid and flexible duct models, the flow field is shown for 18 and 21 generations in Figure 4.The computed flow field is represented by separation streamlines.The obtained flow fields have similar streamlines for both simulated generations.
Figure 4 demonstrates the flow field and pattern 2.5 s after the start of a 5 s inspiration.Figure 4a,b shows the flow field predicted for the 18th generation of rigid and realistic models, respectively, with flow directed from left to right.Similarly, the flow field is characterized by curvilinear streamlines at alveolar openings for the 21st generation in Figure 4c,d.For the rigid model, no bulk convective exchange is observed between the surrounding alveoli and the central part of model; rather, a separation streamline exists at the mouth of the alveoli.In addition, the recirculating flows, which have a velocity several orders of magnitude smaller than the mean lumen velocity, are positioned in the centre of each alveoli segment.Similar flow regimes were determined in previous studies with low Reynolds number flow and rigid boundaries [12,18,19].To satisfy mass conservation, a significant portion of the flow in the central of duct enters the alveoli in the flexible realistic model.Therefore, compared to the rigid model the flow pattern differs, these differences become more tangible with increasing generation number.
A comparison of the recirculating flow velocities in the alveoli for rigid and flexible models shows that the velocity of air flows during the moving boundary model is less than in the non-moving case.In addition, the high velocity region positioned closer to the centreline of the lumen was affected more significantly by the pressure gradient for low generations than high generations.Indeed the

Flow Field and Particle Trajectories
For the rigid and flexible duct models, the flow field is shown for 18 and 21 generations in Figure 4.The computed flow field is represented by separation streamlines.The obtained flow fields have similar streamlines for both simulated generations.
Figure 4 demonstrates the flow field and pattern 2.5 s after the start of a 5 s inspiration.Figure 4a,b shows the flow field predicted for the 18th generation of rigid and realistic models, respectively, with flow directed from left to right.Similarly, the flow field is characterized by curvilinear streamlines at alveolar openings for the 21st generation in Figure 4c,d.For the rigid model, no bulk convective exchange is observed between the surrounding alveoli and the central part of model; rather, a separation streamline exists at the mouth of the alveoli.In addition, the recirculating flows, which have a velocity several orders of magnitude smaller than the mean lumen velocity, are positioned in the centre of each alveoli segment.Similar flow regimes were determined in previous studies with low Reynolds number flow and rigid boundaries [12,18,19].To satisfy mass conservation, a significant portion of the flow in the central of duct enters the alveoli in the flexible realistic model.Therefore, compared to the rigid model the flow pattern differs, these differences become more tangible with increasing generation number.
A comparison of the recirculating flow velocities in the alveoli for rigid and flexible models shows that the velocity of air flows during the moving boundary model is less than in the non-moving case.
In addition, the high velocity region positioned closer to the centreline of the lumen was affected more significantly by the pressure gradient for low generations than high generations.Indeed the streamlines in our model having deformable walls are expected to analyse faraway situations of proximal regions in the acinus accurately.
Figure 5a,b shows the 1 and 5 µm one-way particle trajectories (coloured black), superimposed on the fluid streamlines (coloured with flow velocity magnitude) for the 21st generation and orientation [0, 0, −1] (This means that the gravity acceleration acts on z-direction for the rigid model).It is seen that the 1 µm particle trajectories with small sedimentation velocities are more accurately agree with the flow pathways.In other words, the curvilinear trajectories of 1 µm particles are observed to have better conformity with the fluid streamlines than the 5 µm particles, which diverge from the fluid streamlines.It is clear that 5 µm particles cannot follow the airflow streamlines at the 21st generation.The physical reason of this phenomenon is the large sedimentation velocity of the 5 µm particles relative to the 1 µm particles.Normally inhaled suspension particles due to dilution of concentration have negligible collisions and interactions between particles.Therefore, the use of Euler-Lagrange or Euler-Euler methods can be satisfactory for simulating nano-or micro-particle deposition and transport in the human respiratory systems.If the high velocity gradient, pressure differences and intense particle collisions are challenging conditions, the numerical analysis must be modified to alternative approaches that investigate the fluid-particle interactions.No experimental study is available for investigating one-and two-way coupling of the particle-fluid interaction of micron-scale particulates inhaled into the lung and respiratory system.To the best of our knowledge, only two previous applications of this method in particle dynamic research have been reported [20,21].However these studies investigated particle transport for 3-12 generations.For the efficient visualization of the content, Figure 5b,c compares one-and two-way particle deposition efficiencies for the 21st generations in the rigid models.It can be seen that two-way coupling particles interact with the fluid pathways and have more reliable trajectories than the one-way particles.Finally the outstanding effect of moving walls in alveolar model has been investigated by Figure 5d relative to the similar two-way particle coupling in rigid one (Figure 5c). Figure 6 shows the local particle deposition patterns for 18-22 generations.The concentrated particle deposition observed there demonstrates that direct impaction is the dominant deposition mechanism for Reynolds numbers higher than Re in = 0.13 (G20).Meanwhile, particle-particle interaction effects and local secondary flow effects become significant, leading to more evenly distributed particle deposition patterns, for values lower than Re in = 0.13.
Appl.Sci.2017, 7, 113 9 of 22 streamlines in our model having deformable walls are expected to analyse faraway situations of proximal regions in the acinus accurately.Figure 5a,b shows the 1 and 5 µ m one-way particle trajectories (coloured black), superimposed on the fluid streamlines (coloured with flow velocity magnitude) for the 21st generation and orientation [0, 0, −1] (This means that the gravity acceleration acts on z-direction for the rigid model).It is seen that the 1 µ m particle trajectories with small sedimentation velocities are more accurately agree with the flow pathways.In other words, the curvilinear trajectories of 1 µ m particles are observed to have better conformity with the fluid streamlines than the 5 µ m particles, which diverge from the fluid streamlines.It is clear that 5 µ m particles cannot follow the airflow streamlines at the 21st generation.The physical reason of this phenomenon is the large sedimentation velocity of the 5 µ m particles relative to the 1 µ m particles.Normally inhaled suspension particles due to dilution of concentration have negligible collisions and interactions between particles.Therefore, the use of Euler-Lagrange or Euler-Euler methods can be satisfactory for simulating nano-or micro-particle deposition and transport in the human respiratory systems.If the high velocity gradient, pressure differences and intense particle collisions are challenging conditions, the numerical analysis must be modified to alternative approaches that investigate the fluid-particle interactions.No experimental study is available for investigating one-and two-way coupling of the particle-fluid interaction of micron-scale particulates inhaled into the lung and respiratory system.To the best of our knowledge, only two previous applications of this method in particle dynamic research have been reported [20,21].However these studies investigated particle transport for 3-12 generations.For the efficient visualization of the content, Figure 5b,c compares one-and two-way particle deposition efficiencies for the 21st generations in the rigid models.It can be seen that two-way coupling particles interact with the fluid pathways and have more reliable trajectories than the one-way particles.Finally the outstanding effect of moving walls in alveolar model has been investigated by Figure 5d relative to the similar two-way particle coupling in rigid one (Figure 5c). Figure 6 shows the local particle deposition patterns for 18-22 generations.The concentrated particle deposition observed there demonstrates that direct impaction is the dominant deposition mechanism for Reynolds numbers higher than Rein = 0.13 (G20).Meanwhile, particle-particle interaction effects and local secondary flow effects become significant, leading to more evenly distributed particle deposition patterns, for values lower than Rein = 0.13.The major particle deposition mechanisms in current model lung airways from G18 to G22 can be summarized as:

Comparisons of Deposition Efficiency Predictions Using DPM Approach
We now investigate computationally particle transport by utilizing the one-and two-way DPM methods considering the primary parameters affecting particle transport (each in a separate subsection).

Comparisons of Deposition Efficiency Predictions Using DPM Approach
We now investigate computationally particle transport by utilizing the one-and two-way DPM methods considering the primary parameters affecting particle transport (each in a separate subsection).

Tissue Visco-Elasticity
Many researchers disregard the some important part of alveolar modelling, including both actual particle dynamics and tissue-driven wall motion challenges.Therefore, significant advantages of our 3D developed model are the efficient particle tracking system and the capability of coupling of the FSI technique for modelling the alveolar accurately.Although this approach, which produces the airflow by applying an oscillatory pressure to the tissue, is computationally intensive, it has several desirable features.These features are not available in other studies that used prescribed airflow conditions [22,23].As a result, the wall motion achieved by the elasticity or viscosity of tissue properties can influence the alveolar flow patterns and particle pathways.The relative effect of tissue mechanical properties on particle deposition efficiency can be investigated in two major parts; the elastic and visco-elastic tissues.
Figure 7 shows particle deposition efficiency vs. various tissue elasticity and Figure 8 shows the particle deposition efficiency vs. tissue visco-elasticity values for 1 and 5 µm particle diameter by two-way coupling method on 18-22 generations.For the investigation of the case of elastic tissue, wall deformation is limited due to increasing stiffness (higher E) of tissue, and it causes lower-velocity airflow, which is produced by a negative pressure gradient.Decreasing the airflow velocity causes lower convective forces.Thus, the lower-velocity air flow pulls the particles into the radial alveoli.As a result, the effects of particle inertia and Brownian diffusion are magnified in higher generations by increasing the elasticity of tissue and cause more particles to pass out of streamlines and become deposited on the walls.
In the other case of the visco-elasticity tissue, it is seen that increasing the tissue viscosity does not have any significant effects on the deposition efficiency or particle dynamic transport.The latter is affected by tissue visco-elastic properties, which alter the system's response for changing and loading the convective dynamics of particles.Suitable physiologic values of tissue viscosity (η) are adjusted for providing approximately 0.15 to 1 ratio of tissue stress relaxation time to breathing period [9].According to the results in Figure 6, the variation in time delay of the tissue response in actual conditions, which include gravitational effects, has no significant effect on the particle deposition efficiency and particle dynamics.Under actual conditions, on which the current study is based, the gravitational effects dominate the particle dynamics more than the diffusion effects.In addition, when the size of particles is increased, the particles are driven more strongly by gravity and become less susceptible to diffusion forces [18].Therefore, changing the visco-elasticity or time delay of the tissue response has a milder effect on the obtained results.Furthermore, increasing the generation number causes a lower-velocity airflow produced by the total volumetric flow rate of the generations.Decreasing the airflow velocity causes lower convective forces.By increasing the generation number, the effects of tissue elasticity become more pronounced.Thus, a lower-velocity airflow pulls the particles into the radial alveoli.Note that E = 40,000 dyn/cm 2 and η = 30,000 g/cms are the chosen values for tissue visco-elastic properties in the investigation of other deposition parameters [9].

Tidal Breathing
As shown in Figure 9, it is clear that different breathing treatment periods have noticeable effects on the deposition of smaller particles.This figure shows 1 and 5 µ m two-way particles deposition for different breathing periods.In this part of the study, we investigate how the particle deposition efficiency is affected by breathing time periods.According to a prior study [9], the impactions between particles and walls of alveoli or deposition phenomenon may be obtained when a particle's computed trajectory causes it to pass out of the fluid streamline domain.We have computed the effect of particle diameter and breathing period (λ TB = 5, 7.5, 11.25, 15 s) for 18-22 generations with gravity acting along the z-axis [0, 0, −1].
The results show that small particles, i.e., 1 µ m, have a greater dependence on breathing period than larger particles, i.e., 5 µ m, which experience more sedimentation transport than small particles.In addition, during slower breathing (higher λ TB ), smaller particles are affected more significantly than during rapid breathing (lower λ TB ).This is attributable to the decreasing convection forces during slower breathing (increasing λ TB ); thus, the particles are carried to the deformable wall as a result of diffusion forces.
Thus, an important consequence of these results is noticeable when the smaller particles remain suspended for a long time under study conditions; the particle transport dynamics can be affected significantly by breathing patterns.These trends are shown for 1 and 5 µ m particles in Figure 8.The value λ TB = 5 (s) is selected as the value of tidal breathing period in the investigation of other deposition parameters [9].

Tidal Breathing
As shown in Figure 9, it is clear that different breathing treatment periods have noticeable effects on the deposition of smaller particles.This figure shows 1 and 5 µm two-way particles deposition for different breathing periods.In this part of the study, we investigate how the particle deposition efficiency is affected by breathing time periods.According to a prior study [9], the impactions between particles and walls of alveoli or deposition phenomenon may be obtained when a particle's computed trajectory causes it to pass out of the fluid streamline domain.We have computed the effect of particle diameter and breathing period (λ TB = 5, 7.5, 11.25, 15 s) for 18-22 generations with gravity acting along the z-axis [0, 0, −1].
The results show that small particles, i.e., 1 µm, have a greater dependence on breathing period than larger particles, i.e., 5 µm, which experience more sedimentation transport than small particles.In addition, during slower breathing (higher λ TB ), smaller particles are affected more significantly than during rapid breathing (lower λ TB ).This is attributable to the decreasing convection forces during slower breathing (increasing λ TB ); thus, the particles are carried to the deformable wall as a result of diffusion forces.
Thus, an important consequence of these results is noticeable when the smaller particles remain suspended for a long time under study conditions; the particle transport dynamics can be affected significantly by breathing patterns.These trends are shown for 1 and 5 µm particles in Figure 8.
The value λ TB = 5 (s) is selected as the value of tidal breathing period in the investigation of other deposition parameters [9].

Particle Diameter
Figures 10 and 11 compare 1-5 µ m range one-and two-way particle deposition efficiency for 18-22 generations.In addition, Table 1 is prepared for easier comparing the values for the one and two coupling cases.Deposition efficiencies are displayed as bar diagrams giving total deposition.As expected, there is only a small difference in the deposition predicted for one-and two-way coupled particles.This can be explained by the differentiation among one-way particles coupling, which provides an acceptable approximation of the negligible influence between fluid flow and two-way particles coupling that causes the mutual influence of gas and particles is accounted to fluid momentum equations.Gravitational sedimentation, which causes particles to settle at the bottom of each alveolar cavity, is the main mechanism of deposition for 1-5 µ m diameter particles.As shown in Figure 6, the orientation of the gravitational force is perpendicular to the plan of the model.It is clear that the differing geometries of the alveolar has important role in particle deposition efficiency.The main airflow direction is 150 µ m parallel to the alveolar openings in the developed models.It would be reasonable for the total deposition efficiency to be greater in the model by increasing the distance travelled by particles.Under this condition, particles have more time for gravitational sedimentation to come closer to alveolar walls.Sedimentation is the primary mechanism of deposition in large particles with 5 µ m diameters, because of the very large sedimentation velocities, whereas particle inertia as the Brownian diffusion force dominates the transport of small particles.As expected, it is clear that deposition efficiency of one-and two-way DPM particles is increased by increasing the generation number.Such a phenomenon is due to reducing the Reynolds number and dominant viscous airflow effects in lower (larger number) generations.Moreover, one-way DPM

Particle Diameter
Figures 10 and 11 compare 1-5 µm range one-and two-way particle deposition efficiency for 18-22 generations.In addition, Table 1 is prepared for easier comparing the values for the one and two coupling cases.Deposition efficiencies are displayed as bar diagrams giving total deposition.As expected, there is only a small difference in the deposition predicted for one-and two-way coupled particles.This can be explained by the differentiation among one-way particles coupling, which provides an acceptable approximation of the negligible influence between fluid flow and two-way particles coupling that causes the mutual influence of gas and particles is accounted to fluid momentum equations.Gravitational sedimentation, which causes particles to settle at the bottom of each alveolar cavity, is the main mechanism of deposition for 1-5 µm diameter particles.As shown in Figure 6, the orientation of the gravitational force is perpendicular to the plan of the model.It is clear that the differing geometries of the alveolar has important role in particle deposition efficiency.The main airflow direction is 150 µm parallel to the alveolar openings in the developed models.It would be reasonable for the total deposition efficiency to be greater in the model by increasing the distance travelled by particles.Under this condition, particles have more time for gravitational sedimentation to come closer to alveolar walls.Sedimentation is the primary mechanism of deposition in large particles with 5 µm diameters, because of the very large sedimentation velocities, whereas particle inertia as the Brownian diffusion force dominates the transport of small particles.As expected, it is clear that deposition efficiency of one-and two-way DPM particles is increased by increasing the generation number.Such a phenomenon is due to reducing the Reynolds number and dominant viscous airflow effects in lower (larger number) generations.Moreover, one-way DPM provides higher deposition efficiency than two-way DEM in many generations of the developed model.According to the research of Feng and Kleinstreuer [20], this result is probably due to the particle-fluid interactions, which re-suspend the sediment particles into the airway stream and hence passes from this generation.The two-way DPM method, despite its complexity, provides more extensive particle interaction and transport results than one-way DPM.Therefore, the results of the two-way DPM approach are more reliable than those of the other approach.provides higher deposition efficiency than two-way DEM in many generations of the developed model.According to the research of Feng and Kleinstreuer [20], this result is probably due to the particle-fluid interactions, which re-suspend the sediment particles into the airway stream and hence passes from this generation.The two-way DPM method, despite its complexity, provides more extensive particle interaction and transport results than one-way DPM.Therefore, the results of the two-way DPM approach are more reliable than those of the other approach.

Gravity Orientation
As shown in Figure 11, deposition efficiency, which is a function of generation number and gravity orientation, is investigated with the two-way DPM approach for three orientations: x-axis [−1, 0, 0], y-axis [0, −1, 0] and z-axis [0, 0, −1].Deposition efficiency is defined as the percentage of trapped particles compared to the injected particles to the inlet segment of the model.It is clear that the deposition efficiency increases by increasing the particle size and generation number, regardless of the model orientation.The obtained data demonstrate the significant role of orientation in the amount of deposition and are in agreement with the results of Darquenne and Paiva [18].
The lowest deposition amplitude is computed as gravity acts on the y-axis [0, −1, 0] orientation.Such a phenomenon is due to the same direction of gravity and motion of particles.Indeed, gravity increases the particle velocity and Reynolds number of airflow, which reduces the viscosity forces.The other two orientations, x-axis [−1, 0, 0] and z-axis [0, 0, −1], exhibit approximately similar results because of the radial walls of the alveolar model in this study.The radial walls as shown in Figure 1 are ignored in many previous 1D and 2D studies.However, Darquenne and Paiva [18] show that the radial walls are significant regions of deposition for 0.01-5 µ m particles due to sedimentation.Indeed, the models composed of radial walls modify the efficiency of deposition patterns.

Comparison with Previous Studies
Here the developed numerical model is implemented based on CFD formulations and the FSI technique simultaneously using the ANSYS software, which is useful for developing numerical solutions of more accurate models.
In order to validate the achieved results from the software analysis, we compared them with data and results of similar research assessments.To this end, we compared the results of particle deposition with the results obtained by Darquenne et al. for a duct rigid model [19].This comparison is shown in Table 2, where we specify the same conditions for our model.As can be seen, in some respects, the acceptable agreement is indicated between the data for the two studies.

Gravity Orientation
As shown in Figure 11, deposition efficiency, which is a function of generation number and gravity orientation, is investigated with the two-way DPM approach for three orientations: x-axis [−1, 0, 0], y-axis [0, −1, 0] and z-axis [0, 0, −1].Deposition efficiency is defined as the percentage of trapped particles compared to the injected particles to the inlet segment of the model.It is clear that the deposition efficiency increases by increasing the particle size and generation number, regardless of the model orientation.The obtained data demonstrate the significant role of orientation in the amount of deposition and are in agreement with the results of Darquenne and Paiva [18].
The lowest deposition amplitude is computed as gravity acts on the y-axis [0, −1, 0] orientation.Such a phenomenon is due to the same direction of gravity and motion of particles.Indeed, gravity increases the particle velocity and Reynolds number of airflow, which reduces the viscosity forces.The other two orientations, x-axis [−1, 0, 0] and z-axis [0, 0, −1], exhibit approximately similar results because of the radial walls of the alveolar model in this study.The radial walls as shown in Figure 1 are ignored in many previous 1D and 2D studies.However, Darquenne and Paiva [18] show that the radial walls are significant regions of deposition for 0.01-5 µm particles due to sedimentation.Indeed, the models composed of radial walls modify the efficiency of deposition patterns.

Comparison with Previous Studies
Here the developed numerical model is implemented based on CFD formulations and the FSI technique simultaneously using the ANSYS software, which is useful for developing numerical solutions of more accurate models.
In order to validate the achieved results from the software analysis, we compared them with data and results of similar research assessments.To this end, we compared the results of particle deposition with the results obtained by Darquenne et al. for a duct rigid model [19].This comparison is shown in Table 2, where we specify the same conditions for our model.As can be seen, in some respects, the acceptable agreement is indicated between the data for the two studies.The local deposition distribution changes among the models of different generations with different Re in combinations and values of the velocity ratio (particle settling velocity to mean flow velocity 18µv ) (see Figure 6).There, the total deposition efficiency is observed to increase by increasing the generation number.However, for constant ( V g V ), the total efficiency number decreases by increasing the generation number (see Figure 12).This result is obtained from previous studies [20], and is also observed by Zhang and Kleinstreuer [7].
The local deposition distribution changes among the models of different generations with different Rein combinations and values of the velocity ratio (particle settling velocity to mean flow velocity (  g  ̅ = ρ g  P 2 18µ ̅ ) (see Figure 6).There, the total deposition efficiency is observed to increase by increasing the generation number.However, for constant (  g  ̅ ), the total efficiency number decreases by increasing the generation number (see Figure 12).This result is obtained from previous studies [20], and is also observed by Zhang and Kleinstreuer [7].In order to compare the obtained results with those from previous experimental research, we note that, contrary to the experimental studies, where inhaled particles traverse the tree shape airways and reach the distal lungs, our modelling addresses the transport dynamics of the faction of certain injected particles in an ensemble of randomly.Therefore, we compare our results to higher deposition efficiencies of experimental studies because some particles are deposited in the upper airways and some particles remain entrained but do not reach the distal lungs before being exhaled.Despite this constraint, the published experimental data are consistent with our results, which provide more actual insights for alveolar modelling.During the resolutions, we show that the small In order to compare the obtained results with those from previous experimental research, we note that, contrary to the experimental studies, where inhaled particles traverse the tree shape airways and reach the distal lungs, our modelling addresses the transport dynamics of the faction of certain injected particles in an ensemble of randomly.Therefore, we compare our results to higher deposition efficiencies of experimental studies because some particles are deposited in the upper airways and some particles remain entrained but do not reach the distal lungs before being exhaled.Despite this constraint, the published experimental data are consistent with our results, which provide more actual insights for alveolar modelling.During the resolutions, we show that the small particles having a small characteristic length scale are affected by diffusion-based parameters, i.e., tissue visco-elasticity properties are more intense than for larger particles.This consequence is confirmed by experimental studies [24] that report higher deposition efficiencies.In this study, the inhaled particles are input at the beginning of the lung and traverse more deeply into the alveoli.However, in another study [25] carried out under micro-gravity conditions, the deposition efficiency is higher than the expected values.The physical reason for this phenomenon is the important role of sedimentation effects, since the absence of gravitational forces produce a higher concentration of particles in the deep lungs.
Moreover, the results of breathing periods are consistent with experimental studies showing that slower breathing (increasing breathing periods) may increase the particle deposition efficiency and enhance transport to the distal lungs.In this regard Kim et al. [26] indicate the increase of deposition efficiency by decreasing the volumetric flow rate (slower breathing) for particles 1 < d p < 5 µm in the distal lunges.In addition, Darquenne et al. [24] illustrate the increasing of deposition efficiency with decreasing volumetric flow rate for a whole-lung model by utilizing 0.87-µm particles.We show that small particles are deposited less than the larger particles, which experience rapid deposition in the deep lungs.Besides Kim et al. [26] showed that, when the particle size decreases, the particle deposition distribution shifts toward the distal lungs.

Conclusions
A fully alveolated fluid-structure interaction model has been developed for alveolar breathing mechanics and underlying mechanisms of aerosol deposition in idealised alveolar sacs.We utilized four treatment group breathing periods and twelve visco-elastic properties of lung parenchyma tissue models for investigating the relative importance of the variables affecting deposition such as particle size, gravity orientation, airflow rate or generation number, tissue properties, breathing patterns and particle-fluid interactions.Airflow streamlines, as well as one-and two-way coupling particles deposition efficiencies, have been predicted by the model rather than particle trajectories in the duct model explicitly.The results indicate that modelling the two-way coupling particle system is important because the two-way DPM approach despite its complexity provides more extensive particle interactions and is more reliable than transport results from the one-way DPM approach.The substantial difference between the results of the two approaches is likely due to particle-fluid interactions, which re-suspend the sediment particles in the airway stream and hence pass from the current generation.Some results of the simulations were carried out for different orientations of the model with respect to gravity forces in 18-22 generations of the human lung, and demonstrated the importance of gravity orientation on particle efficiency.Model orientations have important effects on the overall deposition which changed by a factor of approximately 8 between the x-axis [−1, 0, 0] and y-axis [0, −1, 0] orientations for all particle sizes in the 22nd generation simulated.The x-axis [−1, 0, 0] and z-axis [0, 0, −1] orientations have approximately similar results because of the radial walls used in the alveoli model in this study.
The results clearly show that deposition efficiency of one-and two-way DPM particles is increased by increasing the particle diameter and generation number.Such a phenomenon is due to a reduction of the Reynolds number and dominant viscous airflow effect for larger particles as well as lower generations.In addition, the results describe the effects of tissue properties and breathing patterns on particle dynamic transport and deposition efficiency.
Tissue Properties: It is clear that increasing the tissue elasticity (E) causes less wall motion, which lowers flow velocities.Therefore, the value of convention forces that pass the particles out of generations is reduced by decreasing the flow velocities.Thus, the particles are pushed into the radial alveoli by the lower-velocity airflow.As a consequence, the effect of tissue elasticity has been observed explicitly in higher generations because particle inertia and Brownian diffusion is magnified during these generations and causes more particles to pass out of streamlines and become deposited on the walls.It can be inferred that under microgravity (a special condition outside of the human body) the effect of tissue properties can be seen more and more clearly.Therefore, under normal gravity conditions, the propagation depth of the particles is not significantly influenced by tissue mechanical properties because the particles settled to the bottom of the alveolar structure.In addition, during actual conditions, on which the current study is based, gravity effects dominate the particle dynamics with increasing particle size and generation number, and cause rapid particle deposition, which decreases the diffusion forces.Therefore, tissue viscosity has no significant effect on the particle deposition efficiency.
Breathing Patterns: The same results have been obtained for the effect of breathing pattern (different breathing periods) on particle deposition efficiency.The results show that larger particles exhibit little dependence on breathing periods because they experience more sedimentation transport than small particles.Moreover, during slower breathings (higher λ TB ), smaller particles are affected more significantly than during quick breathings (lower λ TB ) because the rapid deposition phenomena and convection forces are weakened during slowed breathing and the particles are carried to the deformable wall as a result of diffusion forces.
The important parameters affecting the particle movement such as particle diameter, tissue visco-elasticity, and tidal breathing period have been optimized with the aid of CFD analysis for a developed model and can be used in pharmacology.By utilizing these results, new drug delivery systems can be improved and modified.In addition, the established analysis system and numerical methods can be applied in magnetic targeting methods for drug delivery to tumour tissue and may also be helpful in designing or modifying drug delivery strategies, producing more effective inhaled pharmaceuticals and treating patients with compromised lung function due to occupational injury or disease.

Limitations of the Model
We know that true breathing patterns in human that include irregular breathing, breath-holding or periodic deep-inspirations, which may not be sinusoidal, and we have not accounted these factors in the simulations.
In actual respiratory system, the streamlines of airflow must enter the alveoli from the lumen by inhalation, and then must exit by exhalation normally.Many studies investigate the convective flow of the lumen and alveolar for demonstrating the effects of expansion and contraction flows.In the model developed here, the particle deposition parameters and flow field properties were only computed during inspiration times over mentioned generations.Therefore, the fate of particles that remain in suspension times or left the model was not computed.In addition, there are anatomical alveolar-alveolar connections that raise the complexity of the gas (and microparticle) delivery.Therefore for simplification of the modelling we ignore some actual anatomical and physiological parameters.
It is expected that by utilizing the realistic model of the human alveolated airways and using similar numerical techniques, more accurate results can be obtained for CFD simulations.

Figure 3 .
Figure 3. Particle deposition percentages for 1-5 µ m range particles as a function of: (a) one-way particle number rate (when the number rate reaches 4000, the deposition percentage becomes stable and the difference between the two steps for increasing of number rate are reasonable and converged to zero); and (b) two-way particles number rate (when the number rate reaches 1500, the deposition percentage becomes stable).

Figure 3 .
Figure 3. Particle deposition percentages for 1-5 µm range particles as a function of: (a) one-way particle number rate (when the number rate reaches 4000, the deposition percentage becomes stable and the difference between the two steps for increasing of number rate are reasonable and converged to zero); and (b) two-way particles number rate (when the number rate reaches 1500, the deposition percentage becomes stable).
(a) The direct inertial impaction causes deposition in the first rows.(b) The 3-D secondary flow effects cause wall deposition.(c) Particle-particle (elastic) contact forces cause wall deposition.
Way Two-Way Diameter One-Way Two-Way Diameter One-Way Two-Way Diameter One-Way Two-Way Diameter One-Way Two-Way

Table 1 .
Comparison between 1 and 5 µ m range particles deposition for generations over 18-22 with gravity acting down in the direction of z-axis [0, 0, −1].

Table 2 .
Comparison of reported and computed results for the same conditions: Percentages of particle deposition in rigid model when gravity force acts along z-axis [0, 0, −1].

Table 2 .
Comparison of reported and computed results for the same conditions: Percentages of particle deposition in rigid model when gravity force acts along z-axis [0, 0, −1].