Fine Particulate Matter Deposition in 3D Out-of-Plane Bifurcation Lung Airway

Received: November 13, 2022
Revised: March 26, 2023
Accepted: April 15, 2023

 Copyright The Author(s). This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are cited.

Cite this article:

Chuang, C.C., Mutuku, J.K., Chueh, C.C., Selvarajoo, A., Chen, W.H. (2023). Fine Particulate Matter Deposition in 3D Out-of-Plane Bifurcation Lung Airway. Aerosol Air Qual. Res. 23, 220392. https://doi.org/10.4209/aaqr.220392

HIGHLIGHTS

• Fine particulate matter (PM_2.5) deposition in lung airways is predicted.
• A 3D out-of-plane bifurcations airways (G8–G12) of Weibel model is developed.
• G8+B_8-9 has the highest deposition fractions as a consequence of high inertial impact.
• The total DFs are inversely proportional to the Stokes number.
• There is an increasing trend for the total DFs as the particle size decreases.

Keywords: Bifurcation airways, Fine particulate matter (PM2.5), Out-of-plane, Computational fluid dynamics (CFD), Aerosol, Particle deposition

1 INTRODUCTION

Fine particulate matter (PM_2.5) containing SARS-CoV-2 would be a perfect carrier. PM_2.5 is the most harmful of all fine particles (Chang et al., 2020; Chen et al., 2021b; Wang et al., 2021). The particulate matter is made up of particle and liquid phases emitted from cooking fumes, road dust, industrial processes, etc. (Lee et al., 2005). The analysis indicated that carbon content and arsenic, selenium, cadmium, antimony, and mercury were the main elements in PM_2.5 found in southern Taiwan (Kuo et al., 2007; Lee et al., 2016) when it came to the biomedicines of human bifurcated bronchi.

Several previous studies have been directed toward connecting the particulate matter (PM) and a passage or airway in the lower respiratory tract that conducts air into the lung (Chen et al., 2021a). Kim and Iglesias (1989) demonstrated one earlier experimental study. They established the correlation between the deposition efficiency and the Stokes number and compared the deposition efficiency associated with different branchial angles. Oldham (2000) completed both experimental and numerical studies, revealing that the numerical results related to the deposition efficiencies of 10 µm particles agreed well with those taken from experiments. In addition, many subsequent numerical studies have focused on pharmaceutical particle delivery or particle deposition in the lung (Sung et al., 2007; Longest et al., 2019; Chalvatzaki et al., 2020). Given that the initial particle deposition arrangements (random-uniform, random-parabolic, and random-random) at the entrance of the inlet, Zhang and Kleinstreuer (2001) showed that the deterministic parabolic distribution had good agreement with the experimental data. Discussing the fate of pharmaceutical drug delivery in the human lung, Sung et al. (2007) found that the particles of porous nanoparticle aggregation might improve the drug delivery for nanoparticles and were likely to avoid mucus clearance effectively. It was shown by the same authors (Zhang and Kleinstreuer, 2002) that the particle deposition becomes stronger with the oscillatory cyclic flow, and the particle deposition distribution of the numerical study is similar to those of the previous experimental studies.

There are two specific human lung diseases: chronic obstructive pulmonary disease (COPD) and asthma. The corollary of the two diseases is a change in the bronchus shape, which presages airflow blockage and breathing-related problems. However, their subtle causes differ between the two diseases. In COPD, the shrinking deformation generally occurs in the middle of the airway, whereas in asthma, the deformation typically takes place on the front of the airway with irregular folds (Chen et al., 2018a; Mutuku and Chen, 2018). Analyzing the flow dynamic in COPD and healthy cases, Mutuku and Chen (2018) found that the contraction of smooth muscles of the COPD bronchus created a physical effect called Dean vorticity. The Dean vortex is induced when the high-velocity fluid encounters a large centrifugal force, ushering in an unstable stratification that makes the high-velocity fluid in the center deflect outwards along the tube bending region. As a result, it forms two counter-rotating roll cells with acute asymmetric flow patterns. Later, utilizing an asthmatic airway and the healthy one, Chen et al. (2018a) showed that the deposition fractions in the asthmatic airway were higher than those in the healthy one due to the significant secondary flow.

A computerized tomography (CT) scan has proven to be a potent tool that combines a series of X-ray images taken from different angles to create the imaging of an in-vitro preparation. The image data are sufficiently manipulated for an in-depth computer-based investigation of the branching pattern of the bronchial tree and statistical geometry analysis. Tracheobronchial airway models have made the simulation results closer to reconstructing the geometry of the bronchial tree of the human lung (Choi et al., 2007; Srivastav et al., 2014; Qi et al., 2018; Dong et al., 2019; Kim et al., 2019). For example, Rahimi-Gorji et al. (2015) used a CT-scan model to create the simulation particle trajectory and deposition pattern. They showed that in the trachea and the subsequent airway for the larger particles (i.e., 5 µm to 10 µm), the deposition fractions significantly increased with an increase in the flow rate, whereas for the smaller particle (i.e., 1 µm or less than 1 µm) the deposition fractions went down when the flow rate increased.

Focusing on the fluid dynamic and particle deposition fractions in tumor tissues in the CT-scan reconstructed lung model, Srivastav et al. (2014) found that the presence of the tumor tissues in the lung ramps up the wall shear stress very rapidly due to the high-velocity gradient. Two optimization methods were used by Chen et al. (2021c) to find out the maximum deposition fractions of the drug particle aerosol in the downstream airway of the asthmatic patient. In doing so, the highest deposition fraction was obtained at the flow rate of 45 L min^–1, the drug dose of 200 µg puff^–1, and the particle diameter of 5 µm with the Taguchi method (Chen et al., 2021c). The use of CT-scan to reconstruct the branching patterns of the human lung remains one of the effectual ways of the mainstream. Using two optimization algorithms, Chen et al. (2021c) analyzed the maximum deposition fractions of the drug particle aerosol in the downstream airway of the asthmatic patient, suggesting that the highest deposition fraction was obtained at the flow rate of 45 L min^–1, the drug dose of 200 µg puff^–1, and the particle diameter of 5 µm via the use of the Taguchi method. In contrast, the maximum deposition fraction was found at the flow rate of 30 L min^–1, the drug dose of 200 µg puff^–1 and the particle diameter of 5 µm through the response surface methodology (RSM).

Past studies have performed numerical calculations for particle deposition fractions using only the geometrical arrangements of bronchioles of three to six consecutive generations, as shown in Table 1. Chen et al. (2021c) solely focused on the flow field and deposition distribution analyses using in-plane configurations of bronchioles, as shown in Fig. 1(a). Later, the flow fields within the bronchi of the generations (G0–G5) following the major bronchi in combination with the out-of-plane configuration were investigated by Pradhan and Guha (2019), as shown in Fig. 1(b). However, they did not consider particle deposition in bronchioles. In the past study, very few studies were carried out on particles deposit within bronchioles of Generations 8 to 12 (G8–G12), especially in the out-of-plane geometry. Therefore, the present study aims to find how the particles deposit within G8–G12 in the out-of-plane geometry. This is useful in showing some early signs of how the particles cause typical human lung diseases within bronchioles. The present article is organized as follows. Section 2 discusses the mathematical modeling for the geometry of bronchioles of G8–G12 of the human lung, the governing equations, and boundary conditions, followed by the results and their discussions in Section 3. The summary and outlook are in the final section.

Fig. 1. Computational domain of G8–G12 of the Weibel lung model for (a) in-plane and (b) out-of-plane configurations.

 
2 PHYSICAL AND MATHEMATICAL DESCRIPTIONS

 
2.1 Bifurcation Geometry of G8–G12

Constructing the geometry of the bronchioles with the different bifurcation angles, physically representing the lung nature, is vital for investigating the effects of acute asthma on the airflow patterns and improving understanding of how the fine particles are transported in response to the flow field inside the lung's bronchioles. For example, the airway tree of some generations would be physically retrieved from CT scans from a single patient (Xu et al., 2020), allowing the assessment of airway wall dimensions that creates the noninvasive investigation of the pathogenesis of airway wall remodeling and the evaluation of new therapeutic interventions. The airway tree can also be made through the Horsfield model, in which Horsfield created an approximate but comprehensive model based on detailed lung castings for human and canine subjects. The Horsfield model of the human lung classified the airway tree into 35 different segment sizes, starting with n = 35, the trachea, and ending with n = 1 (Parker et al., 1971). Nonetheless, using the Horsefield model typically requires specifying the degree of asymmetry at each airway bifurcation, which brings some difficulty in creating the airway tree when parts of the airway tree are considered only. The Weibel model is another model of the respiratory airway tree in which each parent airway, starting with the trachea, splits into two daughter airways. Although the Weibel model is an idealized model (Weibel, 1963), it has proven to offer a good agreement with the experimental data in terms of particle deposition. As a result, instead of using the Horsfield model, for simplicity, the Weibel model is utilized in the present study to construct the airway tree of G8–G12.

 
2.2 Governing Equations for Simulations

The present study assumes that the air flow in the airways is incompressible with a constant viscous coefficient (µ). The flow goes through an isothermal process, meaning the system has no heat exchange. Hofmann et al. (1995) noted that the gravitational effect was pronounced for the fine particles large than 10 µm, and particle deposition fractions increased for particles’ diameters less than 0.1 µm due to the Brownian force.

The gravitational and Brownian effects are neglected in the present study because of the particle diameters between 0.3 and 0.75 µm. The conservation equations of gas-solid two-phase flow are shown as follows.

Continuity equation:

Momentum equation:

where t is the elapsed time, x_i is Cartesian coordinates, U_i is local airflow velocity, ρ is the density of the air, P is pressure, σ_ij is the shear stress, F_pj is the viscous force on the particle. This study aims to analyze the particle motion in response to the fluid flow in the bronchial airways with the momentum interaction between the particle motion and the flow field. The effects of external forces, temperature changes, and magnetic fields are insignificant and hence ignored. The force of the particle motion is expressed as follows.

where m_p is a single particle mass, u_p is the particle velocity, and d_p is the particle diameter. Additionally, there is neither rapture nor deformation on every single particle while tracing with the fluid flow, and the rotation motion is disregarded. In Eq. (3), C_D stands for the drag coefficient based on the combination of Wen and Yu (1966) and Ergun (1952) in the following.

where Re_p is the particle Reynolds number which is expressed as follows.

The concept of void ratio (α_f), representing the scale of the gas phase in the control volume, is applied to distinguish single-phase flow and gas-solid two-phase flow. The ratio is expressed as:

where V_k,p represents the total volume of the k particles in the control volume, and V is the total control volume.

 
2.3 Boundary Conditions

Realistic breathing conditions are crucial in offering accurate information on the bronchial tree's local and regional aerosol depositions. While the actual breathing patterns of each person vary and thus are difficult to define due to their unsteady and irregularity, following Zhang et al. (2002), three typical breathing modes are employed in the present study: resting, light activity, and moderate activity. The three breathing conditions can be approximated through a sinusoidal wave with the respiration rates for inhalation of 15, 30, and 60 L min^–1, respectively, as shown in Fig. 2. Meanwhile, their breathing frequencies are 14, 15.5, and 24.5 cycles per minute, respectively (see Zhang et al. (2002)). The velocity profile at the inlet is position dependent in the following.

Fig. 2. Diagrammatic representations of three breathing conditions for (a) rest, (b) light activity, and (c) moderate activity.

where y is the value of the y-axis coordinate at the inlet, z is the value of the x-axis coordinate at the inlet, R is the radius of the inlet airway. The equation represents the parabolic velocity profile at the inlet, which is considered a fully-developed flow where flow characteristics no longer change with increased distance along the airway. The parabolic profile is, therefore, more suited for the particle trajectory than the uniform profile (Zhang et al., 1997). The average inlet velocity created by Weibel (1963) is applied in this study:

where G is the generation of the lung in Weibel’s geometry, D is the diameter of the entrance bronchial airway, Q_in is the flow rate that enters in the trachea.

 
2.4 Particle Phases

The PM_2.5 concentration of around 55.7 µg m^–3 (Chen et al., 1997) is selected to represent one bottleneck or gridlock where the traffic cannot pass quickly or streets in a city are so full of cars that they cannot move. According to the particle size, the fine particles are typically chosen so that the mode-resolved density of 1.1–2.0 µg m^–3 is achieved, as revealed in Kannosto et al. (2008). As a result, the density of the fine particles of 1.5 µg m^–3 is chosen to ensure the particles' uniform distribution with three particle sizes of 0.3 µm, 0.5 µm, and 0.75 µm considered. The parameters associated with the corresponding numbers of the injected particles under each time step of three volumetric flow rates are listed in Table 2. At the inlet, the particle’s distribution is shown in Fig. 3. A random particle distribution at the inlet is assumed in this study to simulate the reality for the particles entering the lung. The initial velocities of the particles released at the inlet boundary are those coincident with their corresponding flow velocities in place, ensuring that the particles move in a consistent and physical way.

 Fig. 3. Randomly distributed particles at the inlet of G8 with diameters of 0.3, 0.5, and 0.75 µm for (a) rest, (b) light activity, and (c) moderate activity.

It is important to note that the Brownian effect is not taken into account because the considered particles’ diameters herein are always greater than 0.1 µm (Gensdarmes, 2015). The gravitational effect of the particles is also ignored, as the gravitational sedimentation is only dominant in those particles ranging in size from 1 µm to 8 µm (Mutuku et al., 2020). Placed in an environment where the atmospheric pressure of 1 atm is ensured, the whole bronchial tree of G8–G12 of the present study is considered stationary, representing that the shape of the bronchial tree does not change throughout the entire simulation process. Given that the mucus of the airway tract of the lung produced from cells found in mucous glands on the airway walls has fine particles trapped (Inthavong et al., 2010), the fine particle-capturing mechanism is achieved in the present simulations. In other words, upon contact with the wall, particles are typically immediately captured and trapped by the thick mucus, which has adhesive properties. Meanwhile, when the particles travel inside the airway tract without touching the walls of the entire bronchial tree of G8–G12, they leave the bronchial structure in a way that the particles are no longer considered, reducing the consumption of the memory used for the particles left and saving computation resources for those particles.

 
2.5 Numerical Analysis

The commercial software SOLIDWORKS 2018 is used to construct the bronchial airway geometry, which is then applied by the commercial software COMSOL Multiphysics to calculate the internal flow field inside the lung bronchioles of G8–G12. The solid geometry comprises airways tubes and bifurcation sections, as shown in Fig. 4. With a given gas-solid two-phase flow, the two-way coupling accomplished by alternately solving the discrete and continuous phase equations until the solutions in both phases have stopped changing is applied within every time step during the transient simulation. The computational domain inside the bronchial airway geometry consists of several tetrahedral elements. The solver used in the present study is multifrontal massively parallel sparse (MUMS).

Fig. 4. Delineation of the tracheobronchial geometry into (a) airways tubes, (b) bifurcations, and (c) defined sections.

The aerosol dynamics is typically characterized by the deposition fraction (DF) and the Stokes number (Stk). The former is defined as the ratio of the number of aerosols deposited to the total number of those inhaled (Kim and Iglesias, 1989). The latter is defined as the ratio of a characteristic particle time scale to a characteristic flow time scale, quantifying particle inertia's relative importance in a flow. Both physical parameters are expressed mathematically in the following, respectively.

where ρ_d is the particle density, d_p is the particle diameter, U_mean is the average velocity of the flow field, µ is the viscosity coefficient, and D is the airway diameter at G8. Comer et al. (2001) have described a correlation between the deposition fraction and the Stokes number, showing that the particles tend to follow the streamlines smoothly with a small Stoke number (smaller than 0.1), whereas the particles of a larger stoke number (greater than 1) are less sensitive to the flow field, but more likely to be affected by the inertial force, thus causing them to stick to the wall (Lambert et al., 2011).

 
2.6 Grid-independent Test

There are different grid systems with 0.23, 0.83, and 1.67 million (M) elements employed to check the grid independence for average velocities (U_mean) of eight cross-sections during the peak flow of light activity. Fig. 5 depicts the results of the grid-independent tests. According to the tests, the relative error between 0.23 M elements in the mesh and 0.83 M elements in the mesh is 12.1%, and the relative error between 0.83 M elements in the mesh and 1.67 M elements in the mesh is 1.82%. Accordingly, this study selects the second grid system (0.83 M elements) for subsequent simulations and predictions.

Fig. 5. Grid-independent test for three grid systems with the numbers of elements as 230,000, 830,000, and 1,670,000.

 
3 RESULTS AND DISCUSSION

 
3.1 Numerical Validation

The velocity contours in the selected Weibel tracheobronchial airway lung model from G3 to G6 are compared to the work by Lin et al. (2005). The time step, breathing cycle, and flow rate at the geometry's inlet (Generation 3) are 0.02 s, 1.8 s, and 13.8 L min^–1, respectively. This study adopts the same boundary conditions as those applied by Lin et al. (2005). The axial and three cross-sectional velocity contours from Lin et al. (2005) are shown in Figs. 6(a) and 6(b). Meanwhile, the corresponding axial and the cross-sectional velocity contours in the present study are shown in Figs. 6(c) and 6(d). Overall, Fig. 6 depicts that the axial and the cross-sectional velocity contours in the present study are similar to those reported.

Fig. 6. Validation work for the velocity contours for (a) axial velocity distribution and (b) velocity contours for selected cross-sections for the works of Lin et al. (2005), and (c) the axial velocity distribution and (d) velocity contours for selected cross-sections in this study.

The validation of particle deposition patterns in the selected Weibel tracheobronchial lung model from G3 to G5 is performed based on the study by Oldham (2000), which involved an experimental approach and numerical analysis. The particles have a diameter of 10 µm, and the airflow rate at the inlet of the geometry (G3) is 7.5 L min^–1. The deposition patterns from the numerical analysis agree with the investigation by Oldham (2020), as seen in Figs. 7(a) and 7(b). This study has a DF of 76%, which is close to Oldham’s (2020) work of a DF of 80%. Although the particle DFs differ slightly due to the different boundary conditions imposed at the inlet, the degree of accuracy attained herein and the similarity in deposition patterns are close enough for the model to give acceptable accuracy levels. Therefore, the methodology applied here for the fluid flow and the particle tracks are effective for further numerical analysis.

Fig. 7. Validation of the particle tracking system using 10 µm particle deposition patterns in the works of (a) Oldham’s work (Oldham, 2000) and (b) this study.

3.2 Analysis of the Flow Fields

Sinusoidal curves are applied to represent ideal breathing statuses, as presented in Fig. 2. The mean Reynolds number (Re_mean) at Generation 8 (G8) of three ideal breathing statuses are 45.7, 91.4, and 183, as seen in Table 3. A diagram depicting the flow streamlines is applied to describe the fluid motion phenomena in the out-of-plane Weibel lung model. The streamlines in G8–12 during three different breathing statuses of rest, light activity, and moderate activity are depicted in Figs. 8(a–c). Meanwhile, the bifurcation sections between G9 and G10 (i.e., B_9-10) under the three statuses are expanded in Fig. 8(d) for clarity.

Fig. 8. Streamlines of flow in G8–G12 under (a) rest, (c) light activity, and (e) moderate activity, and for B_9-10 under (d) rest, (e) light activity, and (f) moderate exercise breathing status.

The streamlines in G8–12 under the breathing status of rest are depicted in Fig. 8(a), where the bifurcation section between G9 and G10 (i.e., B_9-10) is expanded in Fig. 8(d) for clarity. The streamlines in Fig. 8(d) split smoothly as the cross-sectional area changes at the bifurcation. However, for the light activity shown in Fig. 8(e), streamlines are more chaotic than those of the rest-breathing status, as seen in Fig. 8(d). Furthermore, for the moderate activity shown in Fig. 8(f), the flow field possesses the most significant separation flow motion. Higher flow velocities skew toward the inner wall of G10, and numbers of the low-speed streamlines depart from the primary fluid field and skew toward the outer wall of the bronchial bifurcation (B_9-10) after the airflow collides with the carina region. As the air stream passes through the bifurcation, moderate activity breathing status has the highest minor losses among all the breathing statuses. These results of the streamline patterns have similar findings presented in previous studies when airflow advanced after the bifurcation (Kolanjiyil and Kleinstreuer, 2017; Islam et al., 2019).

3.3 Velocity Contours and Vectors

In the literature, velocity contours and vectors highlighted the dynamics of the fluid motion and its effects on particle deposition (Nowak et al., 2003; Nazridoust and Asgharian, 2008; Chen et al., 2018b). To analyze the airflow phenomena for G8–G12 in the out-of-plane Weibel lung model, some presented cross-sections are defined in Fig. 9, and each section’s velocity contours and secondary flow vectors for the three breathing statuses are presented in Figs. 10–12.

Fig. 9. The location and notation of the out-of-plane G8-G12 Weibel lung model for cross-sections at (a) G9, (b) G10-1, (c) G10-2, (d) G11-1, (e) G11-2, (f) G11-3, (g) G11-4, (h) G12-1, (i) G12-2, (j) G12-3, (k) G12-4, (l) G12-5, (m) G12-6, (n) G12-7, and (o) G12-8.

Fig. 10. Velocity contours and secondary flow vectors under the rest breathing conditions of cross-sections (a) A–A’, (b) B–B’, (c) C–C’, (d) D–D’, (e) E–E’, (f) F–F’, (g) G–G’, (h) H–H’, (i) I–I’, (j) J–J’, (k) K–K’, (l) L–L’, (m) M–M’, (n) N–N’, (o) O–O’, and (p) P–P’.

Fig. 11. Velocity contours and secondary flow vectors under light activity breathing condition for cross-sections (a) A–A’ (b) B–B’, (c) C–C’, (d) D–D’, (e) E–E’, (f) F–F’, (g) G–G’, (h) H–H’, (i) I–I’, (j) J–J’, (k) K–K’, (l) L–L’, (m) M–M’, (n) N–N’, (o) O–O’, and (p) P–P’.

Fig. 12. Velocity contours and secondary flow vectors under moderate exercise breathing condition for cross-sections (a) A–A’, (b) B–B’, (c) C–C’, (d) D–D’, (e) E–E’, (f) F–F’, (g) G–G’, (h) H–H’, (i) I–I’, (j) J–J’, (k) K–K’, (l) L–L’, (m) M–M’, (n) N–N’, (o) cross-section O–O’, and (p) P–P’.

The difference in the total cross-sectional areas between G8 and G9 is 1 mm², while that between G9 and G10 is 1.59 mm². Because of the curvatures at the bifurcations at the end of G8 and G9 (i.e., the B_8-9 and B_9-10), the high-velocity airflow is inclined toward the inner walls (i.e., sides A’, B’, and C’) of cross-sections A–A’, B–B’, and C–C’ shown in Figs. 10(a–c), 11(a–c), and 12(a–c). The secondary flow superposed on the flow field created double vortices with horizontal arrangements at cross-sections A–A’ B–B’, and C–C’. These phenomena can be interpreted as curvatures change in the bifurcations B_8-9 and B_9-10 depicted in Fig. 3(b), which introduces the transverse pressure gradient. The transverse pressure gradient provides the centrifugal force, which results in higher velocities of airflow skewing toward the inner wall of the airways tube. These findings are similar to a previous study where the velocity distributions and particle deposition are studied in gas and particle phase motion in the first and the second bifurcations in a realistic human tracheobronchial airway (Gemci et al., 2007).

The average velocity (U_mean  at the inlets of G11-1, G11-2, G11-3, and G11-4 during the peak flow for light activity are 0.365, 0.337, 0.366, and 0.337 m s^–1, respectively. The velocities in G11-1 and G11-3 are about 8% higher than those in G11-2, and G11-4. Zhang et al. (2002) utilized the in-plane model (G3–G6) to probe the flow pattern and particle deposition pattern in the lung and found that the average velocity (U_average) in the same generation of airway tubes (i.e., G5-1, G5-2, G5-3, and G5-4) relative to the centerline for the in-plane geometry are 1.23, 1.66, 1.42, and 1.04 m s^–1, respectively (Zhang et al., 2002). The velocity in G5-2 is about 35%, 17%, and 63% higher than in G5-1, G5-3, and G5-4, respectively (Zhang et al., 2002). It is evident that the symmetry of the flow field in the out-of-plane configuration is higher than the in-plane geometry. The difference between the average velocity (U_average) and the maximum velocity (U_max) in cross-sections D–D’, E–E’, F–F’, and G–G’ during the peak period of the light activity is 0.135, 0.133, 0.135, and 0.132 m s^–1, respectively. This implies that the initially parabolic flow becomes uniform as the flow advances along the computational domain, as seen in Figs. 11(d–g). This is due to the major and minor losses and loss of inertial forces as the airflow collides with the obstruction at each bifurcation (Xiao et al., 2020). During the rest inhalation status, Dean vortices are not pronounced in the flow fields depicted in Figs. 10(d–g). Owing to the low inertial force of the airflow under the rest, the centrifugal force is inadequate to form the Dean vortices as the airstream traverses the bifurcation (B_10-11). On the contrary, during the light activity and moderate activity, Dean vortices occur in all D–D’, E–E’, F–F’, and G–G’ cross-sections as illustrated in Figs. 11(d–g) and 12(d–g).

During the rest and light activity, the velocity vectors become distorted and move from the periphery area toward the middle zone of airway tubes for cross-sections H–H’, I–I’, and K–K’, as shown in Figs. 10(h), 10(i), and 10(k) and 11(h), 11(i), and 11(k). This is due to the reduction in the influence of the inertial force of the airflow, which is insufficient to sustain the boundary to form Dean vortices, thus resulting in the flow field becoming more distorted. It is noteworthy that Dean vortices appear again in the cross-section J–J’ shown in Figs. 10(j) and 11(j). This indicates that the high-velocity flow fields lie at the airways closer to the centerline of the out-of-plane geometry. The negative centrifugal pressure gradient generated herein is inadequate for generating Dean vortices in the cross-section at G12-3. In contrast, during the moderate activity, velocity vectors demonstrate the occurrence of Dean vortices in all cross-sections at G12 presented in Figs. 12(h–o). Although the major and minor losses increase gradually toward the distal part of the lung, the flow field herein maintains sufficient inertial force to form Dean vortices after the curvatures.

 
3.4 Deposition Patterns for PM_2.5

The PM_2.5 concentration used in this study was 50 µg m^–3, where three PM_2.5 diameters of 0.3 µm, 0.5 µm, and 0.75 µm are applied. These particles are randomly arranged on the X–Y plane and are injected at the entrance of G8 during the respiratory statuses of rest, light activity, and moderate exercise. The goal is to explore PM_2.5 deposition scenarios in G8–G12 of the human lung to effectively point out the deposition hotspots (Nazridoust and Asgharian, 2008; Verbanck et al., 2011; Mutuku et al., 2020) by integrating the airway generations and the bifurcation into a section. The whole geometry is divided into five sections, including G8 + B_8-9, G9 + B_9-10, G10 + B_10-11, G11 + B_11-12, and G12, as shown in Fig. 3(c). The total numbers of PM_2.5 deposited during rest, light activity, and moderate activity are 357, 1049, and 1,351, as seen in Table 4. Table 5 shows the DFs of the five sections. There are nearly no particles deposited in section G12 during all inhalation statuses.

Bifurcation 8–9 (B_8-9) and the periphery of the bronchus at G8 are the particle deposition hotspots, as depicted in Fig. 13. Particles deposit at B_8-9 due to the strong impaction force between the high-velocity gas and particle mixture and B_8-9. The deposition of particles in the periphery of G8 is due to secondary flow fields at G8, as shown in Figs. 10(a), 11(a), and 12(a), where the flow direction diffuses from the central high-velocity area towards the wall. Previous studies (Zhang et al., 2009; Saber and Heydari, 2012; Islam et al., 2017) indicate similar results, where more significant particle deposition occurs in the proximal lung sections compared to the distal ones.

Fig. 13. Deposition patterns for PM_2.5 of diameters 0.3 µm, 0.5 µm, and 0.75 µm in the human bronchial airways (G8–G12) under breathing conditions of (a) rest, (b) light activity, and (c) moderate activity.

During the rest breathing status shown in Fig. 14, the DFs of PM_2.5 in sections G8 + B_8-9, G9 + B_9-10, G10 + B_10-11, and G11 + B_11-12 are 3.49%, 0.26%, 0.17%, and 0.35%, respectively. The low DFs are due to low inertial forces and, therefore, less impaction of particles on the walls. The downstream section of G11 + B_10-11 possesses higher DFs than the upstream sections of G9 + B_9-10 and G10 + B_10-11, despite stronger inertial force. The surface area of section G11 + B_11-12 is responsible for this phenomenon. Even though fewer particles are deposited per unit area, overall DFs exceed those of other sections because of the large surface area compared to sections G9 + B_9-10 and G10 + B_10-11. A similar outcome happens in a study by Feng and Kleinstreuer (2014), where lower lung sections have higher particle deposition than intermediate lung sections.

Fig. 14. The bar chart of the deposition fraction in five defined sections for three selected particle diameters under three inhalation statuses.

For the DFs of PM_2.5 under light activity, the DFs in sections G8 + B_8-9, G9 + B_9-10, G10 + B_10-11, and G11 + B_11-12 are 4.5%, 0.52%, 0.71%, and 0.37%, respectively (Fig. 14). These are increments in DFs of 1.02%, 0.26%, 0.54%, and 0.02% in sections G8 + B_8-9, G9 + B_9-10, G10 + B_10-11, and G11 + B_11-12, respectively compared to those of rest inhalation status. The inertial force increases as the respiratory rate increases; therefore, the impaction of the particles on the walls and the DFs increase. Since B_10-11 has the highest DF (excluding G8 + B_8-9,), it can be inferred that PM_2.5 of the sizes 0.3 µm, 0.5 µm, and 0.75 µm are more affected by the secondary flow fields within the vicinity of B_10-11, including cross-sections D–D', E–E', F–F', and G–G' as seen in Figs. 11(d–g) than those of B_9-10 and B_11-12.

For the moderate activity, the DFs in sections G8 + B_8-9, G9 + B_9-10, G10 + B_10-11, and G11 + B_11-12 are 4.6%, 1.03%, 0.91%, 0.38%, respectively. The DFs increase by 0.1%, 0.51%, 0.2%, and 0.01% in G8 + B_8-9, G9 + B_9-10, G10 + B_10-11, and G11 + B_11-12, respectively, respectively compared to those of light activity (Table 5). Moreover, the DFs tend to decrease toward the distal part of the lung as the intensity of the inertial force decreases in all particle sizes. This is similar to the outcome by Inthavong et al. (2010), which found a significant particle DF at the first bifurcation of the computational domain. Generally, the particles’ DFs decrease as they advance downstream of the lung sections.

The overall trend in the entire geometry indicates that the DFs of PM_2.5 in sections G8 + B_8-9, G9 + B_9-10, G10 + B_10-11, and G11 + B_11-12 increase with the increase in respiratory rates. However, for G11 + B_11-12, as the inhalation status changes from rest to light activity and further to moderate exercise, there are slight increases in DFs of 0.02% and 0.01%, respectively. This behavior indicates that the respiratory status has limited influence on the DFs in the downstream section of the computational domain. An investigation by Inthavong et al. (2010) had similar findings, where particle deposition in a realistic lung model for respiratory rates of 66 and 24 L min^–1 indicated higher variances in DFs for the upstream bifurcations than the downstream ones.

The total DFs and Stokes numbers under different boundary conditions are shown in Table 6. The DFs for the breathing statuses rank as moderate exercise > light activity > rest with respective DFs of 6.92%, 6.13%, and 4.28% for all PM_2.5 (Li et al., 2007; Augusto et al., 2016). As the inhalation status changes from rest to light activity and further to moderate activity, the respective DFs increase by 1.85% and 0.79%, respectively. An increment in inhalation from rest to light activity causes a more significant increase in DF compared to the shift from light activity to moderate exercise (Choi et al., 2007).

Stokes number indicates the ease by which particles get dislodged from their streamlines of flow due to secondary flow motion or their inability to change their flow direction as the shape of the geometry changes. The value of Stokes number at the inlet of G8 ranges from 8.5 × 10^–5 to 2.25 × 10^–3. The correlation between the total deposition fractions and the Stokes number is depicted in Fig. 15. It is evident that the total DFs are inversely proportional to the Stokes number for all inhalation statuses, which contradicts the findings reported in previous studies (Kim and Iglesias, 1989; Li et al., 2007; Augusto et al., 2016). Therefore, it can be inferred that when the particle’s diameter is smaller than 0.5 µm, the Stokes number might not be a good indicator for the total particle DFs in G8–G12 of the out-of-plane lung bifurcation.

Fig. 15. The correlations between the total deposition fractions (DFs) and Stokes number (Stk).

The relationship between the total DFs and PM_2.5 diameters is illustrated in Fig. 16, where the DFs decrease with an increase in particle diameters for all inhalation statuses. This indicates that the effect of flow fields on DFs is inversely proportional to PM_2.5 diameter in the out-of-plane (G8–G12) Weibel airways (U.S. EPA, 2004; Asgharian and Price, 2006; Xi et al., 2008). The total deposition fractions (DFs) observed in this study differ from those in previous studies due to variations in the scales and generations of the lungs being investigated. For example, in a study involving G0–G16, higher DFs were obtained due to the larger computational domain scale (U.S. EPA, 2004; Asgharian and Price, 2006). In another study focused on the upper tracheobronchial region of the lung, G0–G4 showed higher deposition fractions (DFs) due to the turbulent flow and, consequently, increased impaction associated with this region (Xi et al., 2008).

Fig. 16. The correlation between particle diameter and total deposition fraction.

4 CONCLUSIONS

In this study, computational fluid dynamics (CFD) has been applied to predict the deposition of three selected diameters of PM_2.5 in an out-of-plane G8–G12 of the Weibel lung model under three inhalation statuses (i.e., rest, light activity, and moderate activity). The main flow field lies at the airways closer to the centerline of the 3D tracheobronchial lung model. The symmetry of flow in the out-of-plane geometry exceeds that of the in-plane model. The specific deposition fractions (DFs) at different sections increase as the respiratory rate increases, where the highest DFs occurs during moderate activity. Section G8 + B_8-9 has the highest DFs, regardless of respiratory status. This is a consequence of high inertial impact. Except for section G8 + B_8-9, the trend of DFs in the remaining sections has different results during different inhalation statuses. For instance, section G11 + B_11-12 has more significant DFs than the other sections during the rest-breathing status, and section G10 + B_10-11 has more prominent DFs than the other sections during light activity. However, for the moderate exercise mode, the DFs gradually decrease as the flow advances from the inlet to the outlet. The total DFs are inversely proportional to the Stokes number; therefore, it is a poor indicator of the particle DFs for particles whose diameters are less than 0.5 µm. For the correlation between particle sizes and total DFs, there is an increasing trend for the total DFs as the particle size decreases.

NOMENCLATURE

 
Greek Letters

 
Subscript

ACKNOWLEDGMENTS

The authors acknowledge the financial support from National Science and Technology Council, Taiwan, R.O.C., under the grant numbers MOST 111-2221-E-006-102-MY3 and MOST 109-2221-E-006-040-MY3 for this research. This research was also supported in part by Higher Education Sprout Project, Ministry of Education to the Headquarters of University Advancement at National Cheng Kung University(NCKU).