Details of regional particle deposition and airflow structures in a realistic model of human tracheobronchial airways: two-phase flow simulation

Computers in Biology and Medicine 74 (2016) 1–17

Contents lists available at ScienceDirect

Computers in Biology and Medicine journal homepage: www.elsevier.com/locate/cbm

Details of regional particle deposition and airflow structures in a realistic model of human tracheobronchial airways: two-phase flow simulation Mohammad Rahimi-Gorji, Tahereh B. Gorji n, Mofid Gorji-Bandpy Department of Mechanical Engineering, Babol Noshirvani University of Technology, P.O. Box 484, Babol, Iran

art ic l e i nf o

a b s t r a c t

Article history: Received 22 January 2016 Received in revised form 23 April 2016 Accepted 26 April 2016

In the present investigation, detailed two-phase flow modeling of airflow, transport and deposition of micro-particles (1–10 mm) in a realistic tracheobronchial airway geometry based on CT scan images under various breathing conditions (i.e. 10–60 l/min) was considered. Lagrangian particle tracking has been used to investigate the particle deposition patterns in a model comprising mouth up to generation G6 of tracheobronchial airways. The results demonstrated that during all breathing patterns, the maximum velocity change occurred in the narrow throat region (Larynx). Due to implementing a realistic geometry for simulations, many irregularities and bending deflections exist in the airways model. Thereby, at higher inhalation rates, these areas are prone to vortical effects which tend to entrap the inhaled particles. According to the results, deposition fraction has a direct relationship with particle aerodynamic diameter (for dp ¼1–10 mm). Enhancing inhalation flow rate and particle size will largely increase the inertial force and consequently, more particle deposition is evident suggesting that inertial impaction is the dominant deposition mechanism in tracheobronchial airways. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Micro-particle deposition Human tracheobronchial airways CT-scan based images Computational two-phase flow Lagrangian particle tracking

1. Introduction Inhalation and exhalation are the most basic physiological functions required for human body to sustain life. Today, the science of respiration is multidisciplinary; having important contributions from different disciplines of medicine, engineering and basic sciences. There are various processes of human respiratory system such as gas exchange, speech production, particle inhalation and others that can be of interest to researchers from different fields of expertize. The focus of this research is to investigate the air flow behavior and particle deposition in human airways using computational flow modeling. The branching pattern of the airway network is a major factor affecting the deposition of inhaled aerosols. Particle deposition studies in the human airway tree are very significant for characterizing risks of exposure to particulate matter, developing effective drug treatment methods for chronic airway diseases, and calculating dosimetry for occupational hygiene and medicinal treatments [1]. Major challenges need to be overcome to acquire data that is sensitive to lung geometry and function for addressing non-uniformity and asymmetry of particle n

Corresponding author. E-mail addresses: [email protected], [email protected] (M. Rahimi-Gorji), [email protected] (T.B. Gorji). http://dx.doi.org/10.1016/j.compbiomed.2016.04.017 0010-4825/& 2016 Elsevier Ltd. All rights reserved.

deposition in the lungs. These challenges are the representation of realistic airway geometry, the imposition of physiological boundary conditions, and the treatment of turbulence [2]. Accurate 3-D models of the human airway network morphology are of great importance in the study of lung deposition patterns of inhaled aerosols. Knowledge of respiratory aerosol deposition rates and locations is necessary to design effective inhaled medications (i.e., aerosol therapy) that target specific lung regions [3,4]. The earliest studies of airflow in the lung airways include experimental works which presented a few velocity profiles and flow patterns for a double bifurcation model [5,6] followed by numerical studies of velocity field, particle trajectories, and deposition efficiency for 2-D and 3-D single and double bifurcations [7], 2-D steady inspiratory airflow through a double bifurcation model representing three generations of the human central airways [8], reconstructing an average geometrical model of the extra-thoracic airways based on data from computed tomography (CT) scans, magnetic resonance imaging (MRI) scans, and direct observation of living subjects [9] and implementation of a simplified, but more complete model of the upper human airways incorporating a bend between mouth and trachea [10] based on data sets provided by cast of human oral airways [11]. Simulation of micro-particle deposition in human airways have shown substantial differences between a realistic geometry and a simplified one; the realistic geometry provided the best predictions

2

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

p a*, a*1 t g FD K Gprop Yprop

Nomenclature x y z ui ug up IP yþ Rep Rer Rk, Ret dp CT MRI TB CFD DICOM STL LRN FOV CD DE

x coordinate y coordinate z coordinate Mean velocity in tensor notation (m/s) Fluid (air) velocity (m/s) Particle velocity (m/s) Inertial parameter ( Ip = ρdp2Q ) y Dimensionless cell height ( y+ = ν τwall ) ρf f Particle Reynolds number (–) Relative Reynolds number (–) Model constant of k-ω LRN ( Rk = 6, Ret = Particle aerodynamic diameter (mm) Computed tomography Magnetic resonance imaging Tracheobronchial Computational fluid dynamics Digital Imaging and Communication Standard Tessellation Language Low Reynolds number Field of view a a Drag coefficient ( a1 + Re2 + 32 ) r Rer Deposition efficiency (–)

Mean static pressure (pa) Model coefficient, k-ω LRN (–) Time (s) Gravity (m/s2) Drag force (N) Turbulent kinetic energy (–) Transferred product specification (–) Transferred product specification’s loss (–)

Greek symbols

k ) υfω

of regional deposition in comparison with experimental data and hence are needed for accurate evaluation of localized deposition patterns [12–16]. Studies which focused on patient-specific lung geometries have proved to be more pertinent to clinical issues [17,18]. Recently, Rahimi-Gorji et al. [19] demonstrated airflow behavior and particle deposition in three breathing conditions through a realistic human respiratory system (Trachea-G2). In the present investigation, computational fluid dynamics (CFD) was utilized to investigate the effects of several different breathing conditions on the particle (dp¼ 1–10 μm) deposition in a 3-D realistic human airway model extending from mouth to generation G6 of airways. The geometry reconstruction from CT-scan images was performed using 3D-DOCTOR software. The simulated inhalation rates include sedentary (i.e., 10 L/min), light breathing (i.e., 15 L/min), normal breathing (i.e., 30 L/min), semi-heavy breathing (i.e., 45 L/ min) and heavy breathing (i.e., 60 L/min) conditions. For capturing the airflow laminar, transient to turbulent behavior within the airways, the low Reynolds number (LRN) k-ω turbulence model was implemented. One-way coupling was used between the continuous phase (i.e., inhaled air) and discrete phase (i.e., injected particle) to track the particles motion.

2. Model reconstruction and grid generation The airway model includes the oral cavity, pharynx, larynx, trachea, and bifurcating airways representing mouth to G6. The images were obtained from computed tomography (CT) spiral multi-slice scan of a healthy 63 years old, non-smoking male. The slice thicknesses of CT-scan images was set to 0.5 mm and a total of 1023 slices were considered for reconstruction of mouth to G6 using the single-matrix scanner in helical mode, a 500-mm field of view (FOV), 120 kV peak at 50 mA. The CT-scan images – DICOM (Digital Imaging and Communication) files- were imported to 3DDOCTOR image processing software. The slices were segmented and imported with a Standard Tessellation Language (STL) format to CATIA-V5 software. The STL format file was converted there to

Fluid density (kg/m3) Particle density (kg/m3) Kinetic molecular viscosity (μ/ρ) ω Specific dissipation rate (s  1) Kinetic eddy viscosity μT /ρ υT Fluid viscosity μf /ρ υf Βi Model constant of k-ω LRN (–) Kinetic viscosity (Kg m  1 s  1) μ α1, α2, α3 Model constants of the Morsi and Alexander drag law (–) Model constant of k-ω LRN ( σω ¼ σk ¼0.5) (–) σω , σk

ρ ρp υ

( )

(

)

CATPART and then to IGES format. There, The IGES file was imported to ANSYS-Workbench 15 (Fig. 1(a)) to create faces, volumes and extended tubes at inlet and outlets. The morphological specifications of the reconstructed patient-specific geometry are listed in Table 1. Since it is noteworthy to avoid reverse flow at the geometry entryways, they should be extended. The segmented computational domain at this stage is shown in Fig. 1(b). The geometry was meshed by structured tetrahedral mesh using the Meshing module (Fig. 1(c)) and the generated mesh was then read into ANSYS FLUENT 15 for numerical simulation. From the zoomed out parts of Fig. 1(c), it can be seen that the angle of the surface tetrahedral grid is smooth enough to generate a high quality volume mesh. Moreover, the triangular distribution is related to the contour of the configuration. More grid nodes are distributed on the region where the profile has acute changes, which can ensure that the reconstructed geometrical model closely matches the realistic situation. Additionally, to fully capture the turbulent boundary layer profile at the near wall regions, 5 prism layers are added to the geometry surface. Based on the implemented turbulence model, the y þ value for the mesh refinement used was less than 1 (y þ o1). For studying the grid size sensitivity on the flow solution, different grids consisting of approximately 2,410,000 (Mesh 1), 3,350,000 (Mesh 2), 4,250,000 (Mesh) and 5,100,000 (Mesh 4) cells were implemented and axial velocity profiles at two sections (A, B in Fig. 1(a)) were obtained. According to Fig. 2(b), the increase of grid size from 4,250,000 to 5,100,000 cells does not alter the results and the velocity profiles are overlapped. Hence, a mesh with 4,250,000 grid size was adopted for the simulations.

3. Numerical formulation 3.1. Airflow phase The governing equations of the flow phase were considered under steady-state, incompressible and Newtonian fluid assumptions. Although the flow structures and subsequently particle

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

Fig. 1. (a) Anatomical illustration of lung components in the reconstructed model (b) The segmented computational domain (c) 3-D generated volume mesh.

3

4

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

Table 1 The morphological data for the current patient-specific model. Position

Quantity

Inlet diameter Oral cavity

12 (mm) Volume Surface area

7.812 (cm3) 28.138 (cm2) 10.1 (mm)

Length Surface area Mean hydraulic diameter

92.38 (mm) 276 (mm2) 16.57 (mm) 86 (deg)

Left Right

10.31 (mm) 12.2 (mm)

Minimum hydraulic diameter at Larynx Trachea

Branching angle of the main bronchi Main bronchi mean diameter

Fig. 2. (a) Sections of mesh independency study. (b) Grid independency results of the flow study at the two sections (Q ¼15 L/min).

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

5

Table 2 Characteristics of different breathing conditions. Physical state

Sedentary

Light breathing

Normal breathing

Semi-heavy breathing

Heavy breathing

Respiration rate (Q, L/min) Mass flow rate ( ṁ , g/s) Inlet velocity (V, m/s) Reynolds number at inlet (Reinlet) Reynolds number at Larynx (ReLarynx) Particle diameter (μm) Stokes number (St)

10 0. 20417 1.408 994 2175 1–10 0.000416–0.0416

15 0. 30625 2.112 1491 3262 1–10 0.0006245–0.06245

30 0. 40834 4.224 2981 6525 1–10 0.001249–0.1249

45 0. 61251 6.336 4472 9787 1–10 0.0018735–0.18735

60 0. 81668 8.448 5962 13050 1–10 0.002498–0.2498

Fig. 3. (a) Configuration of the Zhao and Lieber model [31] used for validation (b) comparison of velocity profiles.

transport are strongly affected by the transient inhalation/exhalation, the flow fields in steady and unsteady inspiration agree fairly well, except during the deceleration phase and near flow reversal thus the air flow can be assumed steady [13,20]. The associated continuity, Navier–Stokes (Momentum), turbulence kinetic energy and specific dissipation rate equations are given as follows [21]:

Continuity equation:

∂ui = 0, ∂xi Momentum equations:

(1)

6

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

Fig. 4. Comparison of total deposition efficiency for different particle size and inhalation rates as a function of inertial parameter ρdp2Q of present work, Zhang et al. [32], and the empirical curve fit proposed by Stahlhofen et al. [33].

⎡ ⎛ ∂u ∂uj ⎞⎤ ∂ui ∂u ∂ ⎢ 1 ∂p ⎟⎟⎥ + uj i = − + ( υ + υ T )⎜⎜ i + ρ ∂xi ∂t ∂xj ∂xj ⎢⎣ ∂xi ⎠⎥⎦ ⎝ ∂xj

(2)

ui, uj (i, j¼ 1, 2, 3) are the velocity components in x, y and z directions. Turbulent kinetic energy (k) equation:

υ ⎞ ∂k ⎤ ∂k ∂k ∂ ⎡⎛ ⎢ ⎜ υf + T ⎟ ⎥ + Gk − Yk + uj = ∂t σk ⎠ ∂xj ⎥⎦ ∂xj ∂xj ⎢⎣ ⎝

(3)

Pseudo-vorticity (ω) equation:

υ ⎞ ∂ω ⎤ ∂ω ∂ω ∂ ⎡⎛ ⎢ ⎜ υf + T ⎟ ⎥ + Gω − Yω + uj = ⎢ ∂t ∂xj ∂xj ⎣ ⎝ σω ⎠ ∂xj ⎥⎦

3.2. Particle phase

(4)

Also, Gprop and Yprop are transferred product specification and its loss.

k υ T = a* ω

βi = 0.072,

⎛ β − 3 + Ret − Rk ⎞ *⎜ i a* = a∞ ⎟ ⎝ 1 + Ret − Rk ⎠

,

* = 1, a∞

σω = σk = 0.5

Rk = 6,

Ret =

exhalation phase [23]. Due to this fact, the particle deposition during exhalation was not reported herein. For the pressure–velocity coupling, the SIMPLEC algorithm was utilized. The residual values of the governing equations were set to converge when the dimensionless mass residual (total mass residual/mass flow rate) o10  6. After convergence of the steady flow equations, particles were injected in an unsteady state manner. The computations were performed on a PC workstation with 16 GB of RAM and 3 GHz CPUs. Typical run times for steady flow and unsteady particle tracking simulations was approximately 96 h.

(5)

k , υfω (6)

The inflow boundary condition was assumed as mass flow rate, based on the definition of the inlet flow rate according to the five breathing condition (10, 15, 30, 45 and 60 L/min). The inlet flow velocities and corresponded Reynolds numbers for each case are listed in Table 2. The associated boundary conditions were considered as pressure outlet [22] for outlets and no slip condition for walls. In order to satisfy the no-slip boundary condition, the airway walls were considered rigid and the velocity was set to zero on the wall boundaries. Since the wall of respiratory system airway is covered by a mucus layer, the smoothness of the model can resemble the physiological characteristics of the airway. The elimination of this phenomenon would affect the airflow. However, the effect during inhalation is not as noticeable as during the

After the flow field is solved, an analysis of fluid-particle flow was carried out using discrete phase modeling (DPM) by defining the initial position, velocity, size and temperature of individual particles. A stochastic tracking, discrete random walk model with a random eddy life time, was utilized for the turbulent dissipation. DPM determines particle deposition by tracking individual particles as they move through the flow using a Lagrangian approach. Particles were injected at random locations in the inlet plane, relying on an in-house MATLAB code. An injection pattern of a Diskus brand dry powder inhaler (DPI) [24] was implemented as the injection pattern through a text file containing particles with random locations on the inlet surface. Each puff of particle injection lasted for 0.805 s [14] with a time step size of 0.005 s and 400 particles were injected at each time step which resulted in 64,000 total particles injected for each simulation. The trajectory and heat/mass transfer calculations are based on the force balance on the particle and on the convective/radioactive heat and mass transfer from the particle, using the local continuous phase conditions as the particle moves through the flow. The particle trajectory was computed through the equation of the balance of forces acting on that particle. The equation describing the particle velocity, in the Lagrange formulation, in a Cartesian coordinate system has the following form:

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

7

Fig. 5. Velocity magnitude contours and secondary velocity streamlines for different breathing rates in an axial cross-section of (a) oral cavity (Plane 1). (b) Larynx (Plane 2). (c) First bifurcation (Plane 5).

∂up ∂t

= FD(u − up) +

gx(ρp − ρ) ρp

CD = a1 +

, (7)

where FD(u − up) is the drag force per unit particle mass (acceleration due to drag) and

FD =

18μ CDRer , ρp dp2 24

(8)

Here, u is the fluid velocity, up is the particle velocity, μ is the dynamic viscosity of the fluid, ρ is the fluid density, ρp is the density of the particle and dp is the particle aerodynamic diameter. Rer is the relative Reynolds number, which is defined as

Rer =

ρd p u p − u μ

,

(9)

The drag coefficient, CD is obtained using the following formula:

a2 a + 32 , Rer Rer

(10)

where ai are constants that employ to smooth spherical particles over several ranges of Re given by Morsi and Alexander [25]. The turbulent particle tracking was performed by a standard Eddy Interaction Model (EIM), method. A more detailed explanation on this model can be found in [26,,27]. Particles that reach the airway wall were set to trap (which implies that the impacted particles do not bounce; once the distance between particle center and the airway wall becomes less than the particle diameter, the particle adheres to the wall) and entryway faces were set as escape boundaries. One-way coupling is assumed between the air and particle flow fields and the interaction between particles is also neglected since the particle flow is so dilute [28,29]. Injected particles are assumed to be spherical and inert with a density of 1000 kg/m3. Different aerodynamic diameters of particles ranging from 1 μm to 10 μm (1, 3, 5, 7 and 10) are used for the entire simulation. The resulted drug physical properties and concentration

8

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

Fig. 5. (continued)

covers a wide range of dose per puff values of commonly used metered dose inhalers introduced by Terzano [30]. The deposition fraction (DF) of particles was defined as the ratio of the number of deposited particles to the number of those which were injected at the mouth inlet.

4. Results and discussion 4.1. Validation study In order to validate the airflow and particle deposition simulations of the current work, two validation studies were carried out. Fig. 3(a) shows the configuration of the Zhao and Lieber [31]

experimental model that was used for the flow validation of the present study (Re¼ 1036). Comparison of the velocity profiles results and Zhao and Lieber work [31] displayed in Fig. 3(b) reveals that the simulation results are in good agreement with the experimental data. To validate the particle deposition results, deposition efficiencies were computed and compared with Zhang et al. [32] experimental data and the average curve of Stahlhofen et al. [33] based on in-vivo experimental work. In Fig. 4, the total deposition efficiency vs. inertial parameter (IP ¼ ρdp2Q) for

1 μm ≤ dp ≤ 8μm,

30

l min

≤ Q ≤ 90

l min

and

ρp = 0.912 (

g cm3

) is

shown. According to the depicted data, the calculated total deposition efficiency trend is consistent with the results of the experimental data. The IP values of the current work were calculated

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

9

Fig. 5. (continued)

between 151 oIPo 91, 200 (g μm2 s  1) which fits the IP value range of Fig. 4. This implies the validity of the implementation of standard EIM to the present airway geometry. 4.2. Airflow patterns Human tracheobronchial airway is representative of complex tubular systems and hence the local turbulence effects at different flow rates have a significant impact on the airflow structure and particle transport and deposition. Two important parameters that influence the transport and deposition patterns are the local geometry of the tracheobronchial tree and the inhalation patterns. The effect of increasing inhalation flow rate on typical airflow velocity magnitude fields (contours and streamlines) at selected cross-sections; oral cavity (plane 1), Larynx (plane 2), and the first flow divider (plane 5), respectively, is depicted in Fig. 5(a–c). According to the contours, the narrow throat region (Larynx) exhibited the maximum velocity magnitude at all flow rates. As shown in Fig. 5(a), at the upper half of the oral cavity (plane 1), airstreams are deflected and co-rotating vortices are set up and as a result, most of the airflow is passed from the bottom half. The same scenario is evident in Fig. 5(b) at the right half of Larynx (plane 2) where secondary motions are set up due to the geometry deflections and sudden changes experienced in the realistic geometry. As expected, the strength of secondary flow becomes

stronger by enhancing the inhalation flow rate. As illustrated in Fig. 5(c), at the first flow divider (plane 5), the low velocity region is moved toward the outside wall of the first flow divider (main left side lung). The secondary velocity field at this region exhibits distinct and large vortices which are augmented at higher inhalation rates. In Fig. 6, the velocity contours at a sagittal crosssectional view of the airway model (starting from mouth up to first flow divider) are illustrated. By increasing flow rate, the velocity magnitude is continuously augmented. At higher flow rates, it is clearly observed that in the region between the Pharynx (throat) and Larynx, the high velocity region is moved towards the outer wall due to the strong centrifugal forces. Fig. 7 shows the Turbulence Kinetic Energy (TKE) magnitude contour at a sagittal crosssection of airway model for different breathing conditions. It is known that a turbulent laryngeal jet is formed at the larynx where the constriction region causes the airflow to rapidly accelerate and release a jet at the upper region of the trachea that it clearly evidenced in this Figure. The static pressure distribution contours at the airway wall in Fig. 8 indicate that the maximum pressure drop is 165 Pa which occurred in the oral cavity for a flow rate of 60 L/ min. Whereas, in the region from oral cavity up to the first flow divider; the throat region (Larynx) obtained the minimum pressure drop which is attributed to the small cross-sectional area of this particular region.

10

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

Fig. 6. Velocity magnitude contour at a sagittal cross-section of airway model for different breathing conditions.

4.3. Particle motion Particle distribution and transport patterns are useful to illustrate the interaction of airflow structures and particle suspension flow, and hence in order to elucidate physical insight of particle deposition. The 3D views of the local particle transport patterns for particles with aerodynamic diameters dp ¼5 mm at t ¼0.05 s after injection for different inhalation rates are shown in Fig. 9(a). Clearly, as the inhalation rate increases, particle motion is

influenced by flow fluctuations, particles velocity and turbulent dispersion effects are enhanced and more particles can reach the downstream airways. Detailed view of particle transport at t¼0.8 s after injection depicted in Fig. 9(b) illustrates the particle velocity and distribution pattern for normal breathing rate (i.e., 30 L/min). As shown in Fig. 9(b), at higher flow rates, particles motion is highly affected by turbulence dispersion effects. Additionally, due to the combined effects of the generated centrifugal force and the secondary flow induced by the curvature of the main left branch,

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

11

Fig. 7. Turbulence Kinetic Energy (TKE) magnitude contour at a sagittal cross-section of airway model for different breathing conditions.

more particles are directed and tend to be slowed down at the outer wall of this region. It is notable that the contour in Fig. 9 displays the particle velocity magnitude. 4.4. Effect of flow rate and particle size on deposition fractions In order identify the regional deposition of particles; the airway geometry is segmented into 8 zones according to Fig. 3b and local deposition fractions in each zone at various inhalation rates and

particle sizes are calculated and compared in Fig. 10. As inertial impaction is one of the main mechanisms of deposition, it is evident that deposition fraction has direct relationship with particle diameter. However, local deposition fractions vary with inhalation flow rate. At lower flow rates (r15 L/min) where laminar flow dominates, the downstream bifurcations of the main left bronchus experienced the highest deposition rates since the smaller inertial force at lower flow rates, directs more particles toward the lower generations. While, at higher flow rates (Z30 L/min), the onset of

12

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

Fig. 8. Isometric view of pressure distribution across the realistic airway model for different breathing conditions.

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

13

Fig. 9. (a) Locations of particles with dp ¼5 mm at t¼ 0.05 s after injection for different breathing rates. (b) Transport of particles with dp ¼ 5 mm at t ¼0.8 s after injection for normal breathing conditions (flow rate ¼30 L/min).

14

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

Fig. 9. (continued)

turbulence, secondary flow induced convection and curvature of the main left branch will result in enhanced deposition rates in zone 3, 4; directing more particles to the left main branch. Also, the particles depositing in the region of pharynx are more in number than the particles depositing in the region of larynx. Due to the significant increase in particle size; the particle inertia increases largely, so that more particles impact on the laryngeal wall and the particle deposition in the trachea may be resulted from the increase of the particle inertia, the recirculating flow and the turbulent dispersion. To compare the deposition fractions in the left ad right main bronchus, Fig. 11 depicts the regional depositions in left (zones 4, 5 and 6) and right (zones 7 and 8) bronchus at different breathing rates and particle sizes are depicted and compared. This figure

clarifies all descriptions regarding particle deposition in left and right bronchus of human respiratory system. Accordingly, when airflow is laminar (flow rate r15 L/min), deposition fraction in right bronchus is more than left. By increasing the inlet flow rate, more particles tended to deposit in the right main bronchus. As a result, for patients with left bronchus disease, it is recommended to take medication in heavy breathing condition to enhance the percentage of particle deposition in the left bronchus. In Fig. 12, the overall deposition fractions during each breathing rate for various particle sizes are presented. Again, the dominant effect of inertial impaction on deposition; as well as the critical role of turbulent dispersion on particle deposition at higher flow rates is proved.

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

15

Fig. 10. Regional particle deposition fraction in each segmented zone at different particle sizes and inhalation rates.

5. Conclusions The laminar-to-turbulent airflow, transport and deposition of micro-particles were investigated in a realistic model of the human airway from the oral cavity up to generation G6 by two-phase flow simulation. The following conclusions are drawn from the obtained results: (1) Due to the geometry deflections and sudden changes experienced in the realistic geometry, secondary motions are set up which tend to entrap the flowing particles. (2) The flow area constriction at the larynx experienced the maximum velocity change which plays an important role on

the particle transport and deposition at downstream regions. (3) At lower flow rates (r 15 L/min) where laminar flow dominates, downstream bifurcations of the left branch experience the highest deposition rates. (4) By increasing inhalation rate, not only particles motion is highly affected by turbulence dispersion effects, but also due to the centrifugal force induced in the subsequent curvature of the main left branch, more particles are directed towards this region. (5) Based on the results, for patients with left bronchus disease, it is recommended to take medication in heavy breathing condition to enhance the percentage of particle deposition in the left bronchus.

16

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

Fig. 11. Comparison of deposition fraction in left and right bronchus for 1, 5 and 10 μm particles at different inhalation rates.

Acknowledgments The authors wish to express their gratitude to Mr. Hosseinzadeh, medical expert at Parto Mazand Medical Imaging Center, Sari, Iran; for providing us the CT scan images used to generate the airway model.

References

Fig. 12. Variations of overall deposition fractions in the airway model vs. particle size and flow rate.

Conflict of interest The authors declare that there is no potential conflict of interest including any financial, personal or other relationships with other people or organizations within that could inappropriately influence (bias) this work.

[1] A.R. Lambert, C.-L. Lin, E. Mardorf, P. O’shaughnessy, CFD Simulation of contaminant decay for high reynolds flow in a controlled environment, Ann. Occup. Hyg. (2009), mep057. [2] C.-l Lin, M.H. Tawhai, G. McLennan, E. Hoffman, Computational fluid dynamics, Eng. Med. Biol. Mag. IEEE 28 (2009) 25–33. [3] W.H. Finlay, The Mechanics of Inhaled Pharmaceutical Aerosols: An Introduction, Academic Press, London, 2001. [4] T.B. Martonen, C.J. Musante, R.A. Segal, J.D. Schroeter, D. Hwang, M.A. Dolovich, R. Burton, R.M. Spencer, J.S. Fleming, Lung models: strengths and limitations, Respir. Care 45 (2000) 712–736. [5] A.W. Proetz, Air currents in the upper respiratory tract and their clinical importance, Ann. Otol. Rhinol. Laryngol. 60 (1951) 439. [6] R. Schroter, M. Sudlow, Flow patterns in models of the human bronchial airways, Respir. Physiol. 7 (1969) 341–355. [7] J. Lee, J. Goo, M. Chung, Characteristics of inertial deposition in a double bifurcation, J. Aerosol Sci. 27 (1996) 119–138. [8] F. Wilquem, G. Degrez, Numerical modeling of steady inspiratory airflow through a three-generation model of the human central airways, J. Biomech. Eng. 119 (1997) 59–65. [9] K.-W. Stapleton, E. Guentsch, M. Hoskinson, W. Finlay, On the suitability of k–ε turbulence modeling for aerosol deposition in the mouth and throat: a comparison with experiment, J. Aerosol Sci. 31 (2000) 739–749. [10] Z. Zhang, C. Kleinstreuer, C. Kim, Micro-particle transport and deposition in a human oral airway model, J. Aerosol Sci. 33 (2002) 1635–1652. [11] Y.-S. Cheng, Y. Zhou, B.T. Chen, Particle deposition in a cast of human oral airways, Aerosol Sci. Technol. 31 (1999) 286–300. [12] J. Xi, P.W. Longest, Transport and deposition of micro-aerosols in realistic and simplified models of the oral airway, Ann. Biomed. Eng. 35 (2007) 560–581.

M. Rahimi-Gorji et al. / Computers in Biology and Medicine 74 (2016) 1–17

[13] P.W. Longest, M. Hindle, S.D. Choudhuri, J. Xi, Comparison of ambient and spray aerosol deposition in a standard induction port and more realistic mouth–throat geometry, J. Aerosol Sci. 39 (2008) 572–591. [14] K. Inthavong, L.-T. Choi, J. Tu, S. Ding, F. Thien, Micron particle deposition in a tracheobronchial airway model under different breathing conditions, Med. Eng. Phys. 32 (2010) 1198–1212. [15] O. Pourmhran, M. Rahimi-Gorji, M. Gorji-Bandpy, T. Gorji, Simulation of magnetic drug targeting through tracheobronchial airway in presence of an external non-uniform magnetic field using Lagrangian magnetic particle tracking, J. Magn. Magn. Mater. (2015). [16] O. Pourmehran, T.B. Gorji, M. Gorji-Bandpy, Magnetic drug targeting through a realistic model of human tracheobronchial airways using computational fluid and particle dynamics, Biomech. Model. Mechanobiol. (2016) 1–20. [17] C. Fetita, S. Mancini, D. Perchet, F. Prêteux, M. Thiriet, L. Vial, An image-based computational model of oscillatory flow in the proximal part of tracheobronchial trees, Comput. Methods Biomech. Biomed. Eng. 8 (2005) 279–293. [18] F. Farkhadnia, T.B. Gorji, M. Gorji-Bandpy, Airflow, transport and regional deposition of aerosol particles during chronic bronchitis of human central airways, Australas. Phys. Eng. Sci. Med. (2015) 1–16. [19] M. Rahimi-Gorji, O. Pourmehran, M. Gorji-Bandpy, T. Gorji, CFD simulation of airflow behavior and particle transport and deposition in different breathing conditions through the realistic model of human airways, J. Mol. Liq. 209 (2015) 121–133. [20] E. Saber, G. Heydari, Flow patterns and deposition fraction of particles in the range of 0.1–10 μm at trachea and the first third generations under different breathing conditions, Comput. Biol. Med. 42 (2012) 631–638. [21] D.C. Wilcox, Turbulent Modeling for CFD, 2nd ed., DCW Industries, CA, 2002. [22] T. Gemci, V. Ponyavin, Y. Chen, H. Chen, R. Collins, Computational model of

[23]

[24]

[25] [26] [27] [28] [29] [30] [31] [32]

[33]

17

airflow in upper 17 generations of human respiratory tract, J. Biomech. 41 (2008) 2047–2054. H. Bahmanzadeh, O. Abouali, M. Faramarzi, G. Ahmadi, Numerical simulation of airflow and micro-particle deposition in human nasal airway pre-and postvirtual sphenoidotomy surgery, Comput. Biol. Med. 61 (2015) 8–18. M. Broeders, J. Molema, N. Vermue, H.T.M. Folgering, Peak inspiratory flow rate and slope of the inhalation profiles in dry powder inhalers, Eur. Respir. J. 18 (2001) 780–783. S. Morsi, A. Alexander, An investigation of particle trajectories in two-phase flow systems, J. Fluid Mech. 55 (1972) 193–208. D. Graham, P. James, Turbulent dispersion of particles using eddy interaction models, Int. J. Multiph. Flow 22 (1996) 157–175. E.A. Hennick, M. Lightstone, A comparison of stochastic separated flow models for particle dispersion in turbulent flows, Energy Fuels 14 (2000) 95–103. S. Elghobashi, On predicting particle-laden turbulent flows, Appl. Sci. Res. 52 (1994) 309–329. C. Crowe, T. Troutt, J. Chung, Numerical models for two-phase turbulent flows, Ann. Rev. Fluid Mech. 28 (1996) 11–43. C. Terzano, Pressurized metered dose inhalers and add-on devices, Pulm. Pharmacol. Ther. 14 (2001) 351–366. Y. Zhao, B.B. Lieber, Steady inspiratory flow in a model symmetric bifurcation, J. Biomech. Eng. 116 (1994) 488–496. Y. Zhang, W. Finlay, E. Matida, Particle deposition measurements and numerical simulation in a highly idealized mouth–throat, J. Aerosol Sci. 35 (2004) 789–803. W. Stahlhofen, G. Rudolf, A. James, Intercomparison of experimental regional aerosol deposition data, J. Aerosol Med. 2 (1989) 285–308.