Two-phase Flow Dynamics and PM2.5 Deposition in Healthy and Obstructed Human Airways during Inhalation

Received: March 9, 2020
Revised: April 16, 2020
Accepted: April 17, 2020
Download Citation: RIS|BibTeX|https://doi.org/10.4209/aaqr.2020.03.0107  

Cite this article:

Mutuku, J.K., Hou, W.C. and Chen, W.H. (2020). Two-phase Flow Dynamics and PM2.5 Deposition in Healthy and Obstructed Human Airways during Inhalation. Aerosol Air Qual. Res. 20: 1094-1110. https://doi.org/10.4209/aaqr.2020.03.0107

HIGHLIGHTS

• This study numerically investigated two-phase flow in the human airways.
• Inhalation curves for rest, light activity, and moderate exercise were used.
• DFs for different PM_2.5 sizes were 0.12–1.18% in healthy and 0.05–0.49% in COPD.
• Important factors in COPD; jet flows, skewed mass flow rates and dean vortices.
• Depositions at the carina regions could be attributed to inertial impaction.

Keywords: Aerosol; Chronic obstructive pulmonary disease (COPD); Obstructed; Numerical methods; Two-phase flow; Deposition fractions (DFs).

NOMENCLATURE

A             Amplitude of the real breathing curve (Reynolds number)
D             Diameter (mm)
F             Force (N)
G             Generation of Weibel’s airway
C             Carina region of the Weibel’s airway
P             Pressure (Pa)
Q             Flow rate (L s^–1)
R             Radius (mm)
U_i            Airflow velocity (m s^–1)
F_pf           Viscous forces on the fluid molecules.
X_i            Characterizes the Cartesian coordinates,
σ_ij            Shear stress
r              Radial coordinate
Re           Reynolds number
t              Time (s)
U_mean       Average velocity at the airway’s entrance (m s^–1)
U_max        Max velocity at the airway’s entrance (m s^–1)
V             Control volume

Greek Letters

ρ             Density (kg m^–3)
σ             Stress (N m^–2)
µ             Viscosity Coefficient
υ_n            Averaged inlet velocity of parent branch (m s^–1)
θ              Angular coordinate (°)

Subscripts

mean     average value
peak       peak value
I              y-coordinates
j              z coordinates

INTRODUCTION

Scientific evidence shows that obstructive respiratory diseases are mostly caused by exposure of the smooth muscles to irritants such as particulate matter in ambient air which causes contraction of the smooth muscles and subsequent inflammation (Sul et al., 2014). Studies have linked fine particulates (PM_2.5) to a wide range of acute and chronic health defects as well as their exacerbation (Brook et al., 2010; Li et al., 2018; Jiang et al., 2019). Potential pathways for the adverse health effects include; inflammation, oxidative stress, and altered immunity (Ni et al., 2015; Chao et al., 2018; Zhang et al., 2018). A case in point is the mounting evidence suggesting that air pollution from the combustion of fossil fuels is responsible for several adverse health effects in urban areas (Valavanidis et al., 2008; Lelieveld et al., 2015). This information is only available for developed and middle-income countries (Cohen et al., 2005). Model estimations depict a worse situation in developing countries although the true situation remains conclusively unknown (Soriano et al., 2017).

The type of resultant adverse impacts due to inhalation of PM_2.5 is dependent on their chemical composition and mass concentrations (Valavanidis et al., 2008; Brook et al., 2010). The toxicity of particulate matter is inversely proportional to the particle size. As an implication, fine and ultra-fine particulate matters tend to be most toxic among the total suspended solids (Valavanidis et al., 2008). Some studies have shown a steady and autonomous relationship between the concentration of fine particulate matter in the ambient air and the serious health defects caused by air pollution (Li et al., 2018). Further, a statistically significant association has been established between long term exposure to PM_2.5 and reduced life expectancy (Pope III et al., 2009). Specifically, reports on ambient air pollution and mortality show that PM_2.5 is responsible for more than 9% of the global annual mortality rates due to cardiopulmonary disease, cancer of the respiratory system, and acute conditions of the respiratory system in children under the age of 5 years (Cohen et al., 2005). In 2012, the residential areas of more than 80% of the global population was believed to exceed annual WHO AQG for fine particulates of 10 µg m^–3 (Rao et al., 2012).

Environmental protection agencies apply indexes for the enforcement of air quality standards and in efforts to protect the public health from harmful pollutants suspended in the air (Valavanidis et al., 2008). The creation of these indices requires an extensive understanding of the chain of processes involved with air toxics from emission, intake by humans and their physiological effects. Earlier studies have shown a fascinating, but not very conclusive understanding of the possible pathways that associate exposure to aerosols in the ambient air and especially PM_2.5 and mortality caused by cardiovascular disease (Brook et al., 2010). As such, a lot of recent studies are focused on aerosols transportation, transformation and deposition inside the human lungs upon inhalation with ambient air.

The bulk of atmospheric aerosols have a smooth and spherical shape and their density range between 1.1–2 g cm^–3 (Liu et al., 2015). In other studies, PM_2.5 densities in Pittsburg and Beijing were found to be 1.5 ± 0.3 g cm^–3 and 1.41 ± 0.40 g cm^–3 respectively (He et al., 2001; Rees et al., 2004). The variation of velocity with time at the mouth of the nasal cavity can be satisfactorily be represented by a sinewave. Real inhalation curves vary slightly from sinewaves of equivalent amplitude. Previous numerical investigations in the past have used both uniform velocity and parabolic velocity distribution at the inlet, the latter is a better representative of the real situation in the human airways (Lambert et al., 2011; Augusto et al., 2016; Kolanjiyil and Kleinstreuer, 2017).

COPD is characterized by persistent resistance of air-flow in the human airways. This condition is often progressive and the response to inhalation of toxins by COPD patients increases with time (Viegas et al., 1996). Histological studies have characterized this condition with a constriction along the affected airway tubes. The constriction is irregular in human-specific lung models while in idealized geometric models, the constriction is usually axisymmetric (Luo et al., 2007; Chen et al., 2012; Mutuku and Chen, 2018).

The resistance of airflow which is important in the characterization of obstructions in the human lung is represented by the ratio of total pressure drop (Sul et al., 2014). Previous research has shown a direct proportionality between pressure drop and air-flow rate due to amplified secondary flows. The resistance of flow is thus a very important property, especially when dealing with obstructed human airways (Yanai et al., 1992; Kang et al., 2011). There is a significant scale of research on air dynamics, aerosol motion, particle deposition and particle clearance in healthy and obstructed human lungs using numerical methods. Numerical analysis has been used extensively in the past to characterize experimentally inaccessible measurements (Chen, 2001; Chen et al., 2014; Chen et al., 2018b; Zhang et al., 2019). Most of these studies involve multi-phase processes (Chen, 2002; Morel, 2015). In a study by Chen et al. (2012), the deposition fractions (DFs) of pharmaceutical aerosols in obstructed human airways for unsteady inlet conditions exceeded those of a steady one. Studies conducted on obstructed human airways before have concentrated on air flow dynamics, pressure distribution, etiology, and pathogenesis in exacerbation of COPD (Jain et al., 2013; Sul et al., 2014; Ni et al., 2015). However, some important information is still missing in the literature concerning the deposition of PM_2.5 with different particle size distribution in the human airways affected by COPD. This paper is aimed at filling the missing gap by applying numerical analysis tools in the assessment of airflow, aerosol motion, and deposition in healthy and obstructed human airways. The effects of particle size distribution, Reynolds number at the inlet, and geometry on deposition patterns and fractions inside the human lungs are discussed. Exceptional features and essential aspects of results from the numerical analyses are also discussed. Evaluating the deposition patterns and fractions results would give supplementary information needed to understand the exacerbation of COPD as a result of inhalation of PM_2.5.

NUMERICAL FORMULATIONS

Solid Geometry and Mesh of the Bifurcation

The real human lung is very complicated and simplification of the geometry for numerical analysis is the first step. There are two well-known and popular lung geometries for the human lungs; Weibel model which is symmetrical and Horsfield model which is asymmetrical (Weibel, 1965; Horsfield et al., 1971). In a study comparing the particle deposition fractions for the Weibel and Horsfield models found higher deposition fractions in Horsfield's model with respect to the whole lung’s model (Kim and Iglesias, 1989). The Weibel model has higher accuracy and reliability in terms of physiological representation of the human lungs and thus it is more popular among researchers (Weibel, 1963). A comprehensive depiction of the Weibel human lung geometry can be found in other literature (Tsuda et al., 2008; Guha et al., 2016).

The oral cavity and upper regions of the respiratory system are rarely affected by chronic inflammation due to the high efficiency of the lung clearance mechanisms in the upper sections. Therefore, COPD in G5–G8 is one of the few possibilities for a patient affected by airway obstruction and for this reason was adopted for the study. The solid geometries for the healthy and obstructed three-dimensional geometry and mesh grids for G5-G8 of the Weibel’s model were created using SOLIDWORKS and are as depicted in Fig. 1. The solid geometry starts at G5 as the first parent generation which branches into two at a bifurcation angle of 35° with regard to the plane z = 0. The demarcation of the computational domain into grids was based on unstructured tetrahedral elements on the ANSYS FLUENT 19.0 platform. Important parameters in the creation of a physiologically accurate lung model including diameter, length, branching angle, and radius of the curvature are provided in previous works (Chen et al., 2012; Mutuku and Chen, 2018). The carina regions were smoothened after lofting mother and daughter branches using the fillet options. In order to differentiate deposition fractions in the two mechanisms under study here; inertial impaction and complex secondary flow phenomenon, the solid geometries were separated into carina regions and curved tube sections as shown in Fig. 2. The geometry of airways deformed due to COPD was characterized by an axisymmetric constriction at the generation affected by the obstruction. For this study, the narrowest part of the constriction had a 50% reduction in diameter as compared to the unobstructed one, as shown in Fig. 1(b). The extend of obstruction was chosen in a way that the volume reduction and airway surface reduction are consistent with histological studies conducted on obstructed airways (Zhang and Papadakis, 2010; Chen et al., 2012).

Fig. 1. Schematics of (a) block structure for healthy geometry (G5–G8) (b) block structure geometry for G5–G8 wherein G6-1 is affected by COPD (c) mesh for the healthy geometry, and (d) mesh for the airways affected by COPD on G6-2. 

Fig. 2. Labeled branches (G5–G8) and carina regions for (a) healthy geometry, and (b) geometry affected by COPD on G6-1.

Governing Equations

Herein, investigations concerning airflow fields as well as PM_2.5’s deposition patterns and fractions for inhalation during three different statuses were investigated. The air flow was assumed to be incompressible, isothermal, and laminar. The density (ρ) and kinematic viscosity (µ) of the air were assumed to remain constant throughout the inhalation cycle. In addition, the force of gravity on the fine particulates was assumed to be insignificant (Johnson and Jackson, 1987; Zhang et al., 2001). Conservation equations for continuity equation and momentum were applied to govern the gas flow and discrete particle motion.

where, x_i is the coordinate in the cartesian system, t stands for the time elapsed, and U_i and P represent the local velocity of air flow and pressure, respectively. Air density is represented by ρ while σ_ij represents stress. The y and z coordinates are represented by the subscripts i and j . Lastly, F_pf signifies the viscous forces on the fine particulate matter.

The motion of fine particulates in the control volume was governed by Newton's second law of motion which is expressed as:

Assuming the effects of external forces, temperature changes and magnetic field were insignificant and hence ignored, the trajectories of the fine particulate matter can be expressed as:

where an individual PM_2.5’s mass was denoted as m_p, while its diameter and velocity were denoted as d and u_p, respectively. The inertial force exerted by the air-flow to the fine particulate matter was assumed to be the only external force. Further, the fine particulates were assumed to be indestructible and non-distortable under centrifugal forces at the bifurcation as well as the complex secondary flow phenomena in the radial direction of the tube sections.

The equation for drag coefficient (C_d) for fine particulates could be expressed as,

where α₁, α₂, and α₃ are constants obtained through experiments. The PM_2.5’s Reynolds number (Re_N) can be expressed as.

The concept of a fluid volume fraction α_f was used to express the ratio of the volume of the aerosol to that of the air in which it was suspended, and it was expressed as

where the particles, their orders in a sequence, and their volume were denoted by p, k, and V_k,p, respectively. The total volume of the computational domain was denoted as V. For a pure gaseous phase flow α_f = 1.

Boundary Conditions

Velocity Parameters at the Inlet

The flow curve of average air-flow velocity versus time at the inlet was described using a real inhalation curve (Mutuku and Chen, 2018). Real inhalation curves were extrapolated from a study whereby a spirometry test on the velocity variation at the mouth with time had been carried out (Tena et al., 2015). In most of the studies involving aerosol and air motion in the human lungs, the breathing statuses are delineated based on the level of human activity and consequently Reynolds number of flows. The most commonly studied breathing statuses include rest condition (sedentary), light activity, and moderate exercise. For ordinary breathing, the expansion and contraction of the airways are insignificant as compared to the axial airflow velocity, therefore the walls of the airways can be assumed to be rigid (Hughes et al., 1972). The three real inhalation transient curves applied for these investigations are depicted in Fig. 3 which shows the instantaneous and average Reynolds number (Re_mean) variations with inhalation time at the inlet (Zhang and Kleinstreuer, 2002). An equation described by Weibel (1965) for obtaining the average inlet velocity in this case G5 (V_n) was applied in this study

where the generation’s number, diameter, and the airflow volumetric at the inlet are denoted by n, D, and Q_in, respectively.

The velocity profile at the inlet’s circular cross-section was described as a completely formed parabolic flow using a user-defined function (UDF) (Mutuku and Chen, 2018). The airways were smooth and rigid, and a non-slip boundary condition was implemented on the tube’s walls.

Fig. 3. Real inhalation waveforms for (a) rest (b) light activity and (c) moderate exercise.

Injection of Particles and their Motion in the Control Volume

The PM_2.5 particle size distribution applied for this study was obtained from earlier studies (Chen et al., 1997; Lin et al., 2008) and it is shown in Figs. 4(a)–4(c). The fine particulates were assumed to have an average particle density of 1.5 g cm^–3 (Rees et al., 2004). The particles had a deterministic-parabolic distribution at the inlet which was earlier investigated by Zhang and Kleinstreuer (2001) and found to give depositional fractions similar to that of a random distribution. Further information concerning the number of fine particulates inhaled during a single breathing cycle is summarized in Table 1. The ratio of the local concentration of PM_2.5 (C(r)) to the average concentration (C) can be represented by the following equation.

where n is the number of particles in a given ring, for instance r_a ≤ r ≤ r_b. The radius at the inlet and mathematical function for truncating the number of the particles are denoted by R and Int[ ], respectively.

For consistency purposes, the same particle size distribution for a PM_2.5 concentration of 50 µg m^–3 was applied in this study, similar to the one applied by Chen et al. (2018a). Therefore, a constant number of fine particulates were injected at the inlet of the control volume for every time step that is; 6474, 12810 and 25846 for rest, light activity, and moderate exercise, respectively. During the entire numerical calculation, the momentum of the particles was caused by the forces of the surrounding mass of air. The boundary conditions for the discrete phase at the inlet and outlet of the control volume were set as “escape into the control volume” and “escape out of the control volume”, respectively (Chen et al., 2018a). Since the lung’s walls are moist and covered with mucus, the trapping fraction for all the particles that touched the walls was 100%. The gauge pressure was set as zero at all outlets in both the healthy and diseased human airway geometries (Inthavong et al., 2010). PM_2.5 disintegrates only under extreme environmental conditions such as high temperatures and pressures which don’t usually happen in the natural ambient environment or inside the human airways (Zhang et al., 2016). Therefore, in the mathematical formulation, one of the simplifying assumptions is that the particles do not break or undergo any form of deformation during the full inhalation cycle. 

Numerical Analysis

Simulation of two-phase flow in the human airways applies one-way coupling whereby the gas phase affects the solid phase only, but there are no feedback forces from the dispersed phase (Comer et al., 2000; Luo et al., 2007; Rahimi-Gorji et al., 2015). The discrete phase model of CFD is activated by defining several parameters of the dispersed phase, among them their position in the control volume under investigation, velocity, diameter, temperature, mass flow, and time of injection (Rahimi-Gorji et al., 2015). The trajectory calculations use the initial location and parameters of the fluid flow for the calculations.

Conservation equations for both the gas flow and discrete phase motion were solved using ANSYS FLUENT 19.0. Coupled algorithm was adopted because the flow in the region under investigation was laminar and the convergence was determined by pressure-velocity coupling. Previous investigations for instance by Rahimi-Gorji et al. (2015) used the same algorithm. During these investigations, the Navier-Stokes equations and coupled algorithm were solved simultaneously. A three-order MUSCL scheme and Green-Gauss Node were applied as the method of discretization of the fluid’s momentum and base for gradient selection respectively. A standard pressure interpolation scheme and a relaxation coefficient of 0.75 for both velocity and pressure were applied. Before running the numerical simulations, a full breathing cycle was run with the only gaseous flow to get rid of startup effects (Inthavong et al., 2010).

The solutions for the solid phase motion was obtained using a hybrid of Euler and Lagrange approaches provided as the discrete phase model (DPM) in ANSYS FLUENT 19.0 platform. The effects of the inertial force of the continuous phase on the trajectories for the solid phase and their reaction forces were solved for each element in the computational domain. The solution for this couple of conservation equations was attained through solving the equations for the discrete and continuous phase equations alternatively until they both converged. To ensure a smooth breathing curve for all inhalation statuses, a time step of Δt = 0.035 s was selected for solving the air-flow fields and paths of the suspended particles. For a single time-step, the solutions of the flow fields converged when the continuity residual had a relative error of less than 10^–5.

Dimensionless Numbers

The inlet’s Reynolds number (Re_d) is related to airflow velocity at the geometry’s inlet as shown below (Verbanck et al., 2011)

where ρ is the fluid’s density (=1.225 kg m^–3), μ represents the viscosity coefficient of air (1.789 ´ 10^–5 kg m^–1 s^–1), and D is the diameter at the inlet of the parent tube (G5) which is 3.5 ´ 10^–3 m. The flow velocity at the entrance of G0 (i.e., the mouth’s entrance) is uniform. But as the flow proceeds deeper into the lower generations, the velocity profile becomes parabolic distribution due to the drag force on walls of the airways. The velocity distribution at the inlet is parabolic, and the maximum flow velocity (U_MAX) at the middle of the cross-section is twice the U_mean. Since the fluid’s density and the viscous coefficient are fixed, the larger the U_mean and U_MAX, the greater the Reynolds number, hence the greater the inertial force of the suspended particles.

The formula for deposition fraction (DF) represents the number of particles deposited as a percentage of the total number of particles entering a particular section of the control volume (Kim and Iglesias, 1989), and is defined as:

For a fixed number of particles in the tube, the larger the DF, the higher the number of particles deposited.

Model Validation and Grid Independence

The axial flow velocity diagrams, velocity contours, and secondary flow vectors for selected cross-sections as well as deposition patterns during inhalation were validated whereby the air flow rate for the inlet specified case was 250 mL s^–1. Schematics for the validation of axial air velocity contours as well as secondary vectors are shown in Fig. 5(a). From the cross-sectional velocity and secondary vector distribution diagrams, it was observed that there were four vortices generated in all the three cross-sections, which was also consistent with the phenomenon observed in an earlier study (Nazridoust and Asgharian, 2008). The overall trends for deposition from this study are consistent with the literature, as shown in Fig. 5(b), and the particles’ deposition patterns were significantly influenced by the different boundary conditions imposed on the breathing cycles. The deposition fractions for 0.1, 1 and 10 µm particles were 7.4, 0.11, and 0.029%, respectively. This closely agreed with earlier findings (Nazridoust and Asgharian, 2008), thereby confirming the accuracy and confidence of the developed model for the prediction of PM_2.5 motion and deposition patterns as well as fractions. Through the above verification, it was clear that the results from the coupled algorithm can be applied to predict air and particle motion in the human airways (Chen et al., 2018a). Therefore, the coupled algorithm was used to simulate the instantaneous two-phase flow in the multiple bifurcation respiratory tracts.

Three different grids with 0.44 million, 1.37 million and 1.87 million elements were applied to check the grid’s independence for the deposition fraction of the particles in a COPD geometry during inhalation for light activity as shown in Fig. 6. With 12806 particles inhaled during the light activity which was characterized by 57 time-steps of 0.035 s each, it was discovered that 1.37 million elements in the mesh were enough to capture the particle deposition patterns and fractions within a reasonable time frame and without significantly compromising on accuracy and reliability of the numerical analysis. Similarly, 1.37 million elements were adopted for simulating the healthy case. 

Fig. 5. Validation for (a) axial velocity distribution and (b)Validation for particle deposition with three different diameter sizes 10 µm, 1 µm and 0.1 µm using results from an investigation by Nazridoust and Asgharian (2008). 

Fig. 6. Grid independence test for the COPD mesh.

RESULTS AND DISCUSSION

Fluid Flow Dynamics

Axial Airflow Fields

Airflow fields are very key in determining the deposition patterns of particulate matter in human airways. Several parameters can be used to explain fluid flow in the human airways such as velocity distribution, stagnation zones, recirculation zones, jet flows, and shear stress on the walls (Zhang et al., 2002; Kleinstreuer and Zhang, 2010; Chen et al., 2012). The airflow velocity distribution affects pressure distribution which was an important parameter in determining hotspots for aerosol deposition (Luo et al., 2007; Kang et al., 2011). The velocity distribution inside the control volume was responsible for the formation of stagnation and recirculation zones as well as jet flow phenomena in the constriction as shown in Fig. 7 (Chen et al., 2012; Chen et al., 2018a). The jets were extracted at the cross-section in the middle of each generation under investigation. These flow phenomena, as a result of velocity distribution, greatly affected particle deposition patterns in both healthy and deformed human airways. Jets tend to form in obstructed sections according to the results from earlier analysis (Kleinstreuer and Zhang, 2010; Sul et al., 2014). Boundary layer separation in the regions near the constrictions leads to the recirculation phenomenon as it can be seen in the vicinity of the obstruction in Fig 7(b) (Sul et al., 2014). Further complex secondary flow phenomena can lead to the formation of Dean vortices. Prior to investigating the two-phase flow, detailed investigations and discussions for the pure gas-phase flow for the geometries discussed here were performed (Mutuku and Chen, 2018). 

Fig. 7. Stream lines of flow and jet flow phenomenon for (a) healthy human airways and (b) COPD geometry during (i) rest, (ii) light activity and, (iii) moderate exercise.

Wall Shear Stress

Wall shear stress is an important factor in the immune response mechanism of the lungs through the release of ATP, which increases the intracellular concentration of calcium for cell engulfment and anti-inflammatory response (Gronski et al., 2009). Fig. 8 depicts that the wall shear stress was maximum in the vicinity of the carina regions as well as around the constriction in the geometry affected by COPD. Jet flow phenomenon in the geometry affected by COPD induced high-velocity gradients in the radial direction. Since wall shear stress is proportional to the velocity gradients perpendicular to the axial airflow direction, the location of the maximum shear stress adhered to the regions where high axial velocities were greatly skewed. Low shear stresses (below 0.3 Pa) have a positive impact on the barrier function of the epithelial lining of the airway walls (Sul et al., 2014). On the other hand, highly exaggerated shear stresses can damage the epithelial cells.

During rest, the healthy generation had shear stress exceeding 0.3 Pa at carina 5 (C5) only, as compared to the COPD geometry whereby exaggerated shear stresses were present at C5, the obstructed generation (G6-1), as well as the carina 6-2 (C6-2). This was because high resistance of flow at the obstructed section led to a higher flow rate in the unobstructed generation. Therefore, high shear stress was generated once the amplified velocity of flow impacted on the carina region. During the light activity, all carina regions except C7-3 and C7-4 had high shear stresses. Additionally, the wall area in the vicinity of the constriction with the highest shear stresses is larger than that observed in the same geometry during rest. The regions of airway walls affected by high shear stress during moderate exercise corresponded to those of light activity but with a larger effective wall area for the hotspots, as shown in Fig. 8. Overall, the geometry affected by COPD had more regions of high shear stress compared to the healthy human airway, indicating a high likelihood of epithelial damage in the obstructed human airways. 

Fig. 8. Wall shear stresses for (a) healthy human airways and (b) COPD geometry during (i) rest, (ii) light activity and, (iii) moderate exercise.

Deposition Patterns for PM_2.5

Three real and transient inhalation curves with an average Reynolds number of 164.3, 362.4, and 606.4 were applied for these investigations as shown in Fig. 3. A PM_2.5 concentration of 50 µg m^–3 was applied where particles with four different diameters including 0.075 µm, 0.15 µm, 0.3 µm, and 0.6 µm as presented in Fig. 4 were injected at the inlet of G5. The PM_2.5 particles exceeding 1 µm were excluded from these investigations because they are easily deposited in the oral, nasal and upper sections of the human airways. In this simulation, G5 is the inlet of the control volume and only particles with less than 1 µm in diameter reach this section of the respiratory system. The results presented herein depict the deposition fractions and patterns obtained from tracking a total of 350031, 692596 and 833553 PM_2.5 particles during rest, light activity, and moderate exercise inhalation statuses, respectively. 

Fig. 4. Distribution of PM_2.5 injected at the inlet for each time step during (a) rest and (b) light activity and (c) moderate exercise.

There are seven main deposition mechanisms for aerosols in the human airways including (1) turbulent mixing; (2) inertial impaction; (3) impaction due to complex secondary flows; (4) centrifugal forces; (5) gravitational sedimentation, (6) Brownian motion; and (7) electrostatic precipitation (Finlay and Martin, 2008; Darquenne, 2012). The dominance of any of these mechanisms of deposition varies as the particles advance from the oral cavity, through the upper section of the human lungs, central, and the lower sections, and later into the alveolar region. The effect of gravitational sedimentation in the overall deposition is usually more pronounced for aerosols exceeding 10 µm (Hofmann et al., 1995). The present study focused on particles with diameters smaller than 2.5 µm. Therefore, the effect of gravity on particles was ignored. Additionally, diffusional forces are more dominant in the alveolar section of the human lungs only and therefore was ignored for this study. In these investigations, the impaction of the PM_2.5 on the airway walls due to inertial forces, centrifugal forces, and complex secondary flow phenomena were investigated. The most important obstacles for inertial impaction in this study are the carina regions. Therefore, to study the hotspots for deposition, both the healthy and obstructed geometries were segmented into; carina region and other sections of the control volumes herein referred to as tube walls as depicted in Fig. 2.

Several factors influencing the deposition fractions and patterns of inhaled particulate matter include (1) geometry of the tracheobronchial model; (2) inhalation status; and (3) chemical and physical properties of the particles (Daigle et al., 2003; Zhang and Papadakis, 2010; Verbanck et al., 2011; Li et al., 2018). The chemical properties of PM_2.5 were not considered for this study; therefore, the particles were assumed to be inert. Since they had a constant density, the PM_2.5 diameters can be directly linked to their Stokes’ numbers. Similar to previous investigations, the results showed that deposition fractions for aerosols are dependent on their Stokes numbers since inertial impaction is the main deposition mechanism in the human airways (Cheng et al., 1999; Zhang et al., 2002). As it can be seen in Table 2, deposition fractions in the carina regions dominated DFs in the whole geometry. The fine particulates under investigation here were assumed to have a spherical shape and consequently, their Stokes number increased with an increase in diameter. 

Deposition on the Carina Regions

Aerosols deposit on carina regions due to abrupt geometrical changes. Therefore, inertial impaction was the main mechanism credited for deposition (Huang and Zhang, 2011). Herein, the stokes number was the ratio of PM_2.5 stopping distance to characteristic distance from the carina regions (Wessel and Righi, 1988). As such fine particulates with higher stokes number tend to deviate from the streamlines of airflow at the carina region, while low stokes number particulates adhere to the streamlines of flow.

As can be seen in the graphic representation of deposition patterns in Figs. 9 and 10, the most important hot spot for PM_2.5 deposition from this study was carina 5 (C5). It is the first carina region in this control volume, and it received the highest DF of fine particulates, as seen in Table. 2. The property of deviation of PM_2.5 from the streamlines of flow in the two-phase flow under investigation is best described using the stokes number which is expressed as;

where St is the Stokes number, d_p and ρ_p are the diameter and density of the PM_2.5, respectively; u and μ are properties of the gas phase that are the average velocity and dynamic viscosity, respectively, and  is the diameter of the airway. Overall, deposition by inertial impaction increases with an increase in the Stokes number at the inlet of the geometry. For all the carina regions (C5 to C7-4) deposition fractions for the healthy geometry exceeded that of COPD geometry for 12 out of 21 instances. Overall, deposition fractions were highest on carina 5 and lowest for carina 7-1 and 7-4 for both geometries. Important factors for the deposition of the PM_2.5 on the carina regions include (1) the position in the control volume, (2) airflow velocity, and (3) particle distribution at the inlet which was parabolic deterministic in this case (Zhang and Kleinstreuer, 2001). Distal carina regions had lower deposition fractions as compared to the proximal ones.

Since the effective cross-sectional area increases as air advances from the inlet (G5) to the outlet (G8), the velocity decreases. This consequently led to a drop in the inertial impaction forces as the particles progress from the inlet to the outlet. For carina 5 and 7-2, particle deposition in the healthy human airway exceeded that in the COPD geometry for all inhalation statuses. In carina 6-1 ad 6-2, particle deposition in the COPD geometry exceeded that in the healthy geometry during all inhalation statuses. In carina 6-2, the deposition fraction of PM_2.5 for the COPD geometry exceeded that deposited the healthy geometry by 22.33, 231.36, and 248.52% for rest, light activity, and moderate exercise respectively. In carina 6-1 which succeeded the obstructed generation, total deposition fractions in the COPD geometry exceeded the ones in the healthy generation for all inhalation statuses by 137.86, 77.05, and 47.45% for rest, light activity, and moderate exercise, respectively. This can be attributed to jet flow phenomena at the obstructed generation. High axial velocities associated with the jet flow led to amplified DFs at the succeeding carina region (C6-1). Due to the high resistance of flow in the obstructed generation (G6-1), the axial velocity of flow in G6-2 was higher compared to the corresponding generation in a healthy geometry. Consequently, the deposition in C6-2 and C6-1 in COPD human geometry exceeded those in a healthy geometry for all inhalation statuses.

Overall PM_2.5 deposition fractions at all the carina regions of the healthy human airway geometry exceeded that of COPD by 9.28, 18.20, and 13.05% for rest, light activity, and moderate exercise respectively. Fig. 11 depicts a scatter of the deposition fractions against the particles’ Stokes numbers at the carina regions, tube walls, and overall deposition at the end of the full inhalation cycle for a COPD and healthy geometry. From the scatter in Fig. 11(a), it is clear that deposition in the carina region was more strongly related to the Stokes number in a healthy geometry than in a COPD geometry. Overall, the deposition of PM_2.5 in the carina regions was more important as it was responsible for 64–88% of the overall deposition fraction. 

Deposition along the Tube Walls

Even though deposition along the tube walls was mostly attributed to complex secondary flow phenomenon and centrifugal forces, inertial impaction could not be fully ruled out. Indeed, deposition at the inner walls of the bifurcation can be attributed to inertial impaction. There was a jet flow phenomenon in the generation affected by COPD (G6-1), and hence the deposition fractions at this generation were very low as compared to its corresponding generation in the healthy geometry. The presence of an obstruction in G6-1 leads to amplified PM_2.5 deposition in G7-2, as shown in Fig. 9. As a result, there was a significant asymmetry in the deposition patterns in G7-2 and G7-3 unlike in the healthy human airway geometry where deposition in G7-3 slightly exceeded that of G7-2. The second most important obstacle in a human airway bifurcation besides the carina regions were the curved sections of the airway tubes. Centrifugal forces played a key role in the axial velocity distribution at the curvatures. As such, there was a slightly amplified deposition of fine particulates on the inner tube walls immediately past the carina regions.

Fig. 9. PM_2.5 deposition in healthy human airway bifurcation (G5– G8) at the (i) carina region, (ii) tube walls and (iii) whole geometry during (a) rest, (b) light activity and, (c) moderate exercise.

Overall Deposition

PM_2.5 deposition fractions in the healthy human geometry exceeded deposition in a COPD human airway geometry during all inhalation statuses except during moderate exercise. The total deposition fractions of PM_2.5 for the healthy geometry were 1.22, 2.50 and, 1.86%; while that of the COPD geometry were 0.51, 0.91and, 2.11 for rest, light activity and moderate exercise respectively. The relative disparity in the deposition fractions of PM_2.5 between the healthy and COPD case for PM_2.5 was 59.19 and 33.30% during rest and light activity respectively. On the other hand, during moderate exercise, deposition fractions in the COPD geometry exceeded that in the healthy geometry by 0.37%. From the graphical representation of deposition fractions in Fig. 12, it can be concluded that the dominance of deposition mechanisms under study namely; inertial impaction, centrifugal forces and complex secondary flow phenomena varied differently for a particular diameter size at different average Reynolds numbers. The particle distribution at the inlet also played a role as particles near the airway walls tend to deposit more easily at the slightest deviation from the streamlines in the axial direction of flow. Further, particles at the rings that are closest to the middle of the cross-section deposited easily upon impaction at the carina region. For a healthy human airway geometry, deposition fractions increased as the Reynolds number increased from rest to light activity. However, as the Reynolds number of flows increased from light activity to moderate exercise, the DFs for the 3 and 6 µm particle diameters decreased slightly. This can be explained by an exaggerated jet flow phenomenon throughout the bifurcation which purged the particulate matter to the exit of the control volume at G8. The stokes numbers of the PM_2.5 particles at the inlet of G5 ranged from 2.012 × 0^–8 to 5.55 × 10^–5. A graphical representation of the overall deposition patterns of the particles is shown in Fig. 10, whereby deposition fractions at the carina region together with overall PM_2.5 deposition fractions could be closely linked to the stokes numbers of the particulate matters, as shown in Figs 10(a) and 10(c) respectively. This is however not the case for the deposition fraction in the tube section as shown in Fig. 10(b) which are mostly due to complex secondary flows and whose presence is independent of the particles’ Stokes numbers. 

Fig. 10. PM_2.5 deposition in an obstructed human airway bifurcation (G5– G8) at the (i) carina region, (ii) tube walls and (iii) whole geometry during (a) rest, (b) light activity and, (c) moderate exercise. 

Fig. 11. A scatter for the deposition fractions vs Stokes number for the (a) carina regions, (b) tube walls and, (c) whole geometry. 

Fig. 12. Overall deposition fractions in (a) healthy human airway geometry and (b) human airway geometry affected by COPD.

CONCLUSIONS

CFD was applied to investigate airflow and localized deposition of fine particulate matters in healthy and obstructed human airways. Herein, complex physical phenomena were broken down and derived from otherwise experimentally inaccessible regions. Wall shear stresses in COPD geometry were found to exceed those in the healthy human airway geometry. The carina regions were responsible for 64–88% of overall PM_2.5 deposition. The most important hotspot for the deposition fine particulate matters in this study was found to be carina 5. Jet flow at the constriction affected by COPD led to amplified DFs at the succeeding carina region. The presence of an obstruction in G6-1 leads to a more amplified PM_2.5 deposition in G7-2 for the COPD as compared to the corresponding generation in the healthy human airway geometry. The intensity of the stagnation and recirculation zones increased with an increase in Reynolds number from 164.3–606.4. Overall PM_2.5 deposition in the carina regions as well as the curved sections was due to inertial impaction, complex secondary flow phenomena and centrifugal forces. PM_2.5 deposition in the human airway affected by COPD was skewed. Overall, the DFs of PM_2.5 increased with an increase in Stokes number. Total deposition fraction of PM_2.5 in a COPD geometry exceeded that of a healthy geometry during the moderate exercise inhalation status only. However, the direct correlation between DFs and stokes number was lacking in a few instances due to very low ranges of Stokes number.

ACKNOWLEDGMENTS

The authors acknowledge financial support from the Ministry of Science and Technology Taiwan, ROC, under the grant numbers MOST 108-2622-E-006-017-CC1 and 109-3116-F-006-016-CC1 for this research.