Numerical Simulation of Transport and Deposition of Dust Particles in Human Tracheobronchial Airways

Dust is a common pollutant of the air we breath. If dust particles are inhaled and deposited in human airways, they can cause a variety of respiratory disorders. The inhaled dust particles motion in human airways goes along with the airflow. The transport process can be considered as a two-phase flow of a gas phase and a particle phase. In this study, we investigated the airflow and dust particles transport and deposition in human tracheobronchial airways using computational fluid dynamics (CFD) techniques. A steady simulation was performed in asymmetric tracheobronchial airway mode consisting of 19 outlets to observe the characteristics of airflow fields. The discrete phase model (DPM) was employed to predict the particle trajectories and deposition in the airway model. Deposition resulted from inertial impaction and gravitational sedimentation was considered. In the simulation, the airflow characteristics differences in the right and left bronchial trees were observed. The influence of secondary flow on dust particles motion was great. More dust particles were deposited in the right bronchial tree than in the left. The deposition fraction of dust particles in human tracheobronchial airways was high. This study can provide awareness on deposition of dust particles passing beyond the larynx and enhance prevention of their entry into the respiratory system. It can also contribute a convenient way on the location of deposition of particles of a given type in human respiratory tract to be used for respiratory disease preventions.


Introduction
Dust is a typical air contaminant generated by many different sources and human activities. Some of the sources of dust are crushing and grinding of rock and mineral products, abrasion of surfaces, farming, non-paved roads and natural erosion of soil [1]. Dust particles suspended in air are likely to be inhaled during breathing. If dust particles are inhaled and deposited in human airways, they can cause a variety of respiratory disorders [2]. So, it is important to understand the transport and deposition of dust particles in the human respiratory tract to treat the respiratory diseases caused by dust pollution from different sources and activities.
Numerous CFD simulations of gas-particle flow have been performed in human airway models using Eulerian-Lagrangian approach [3][4][5][6][7] to investigate the transport and deposition of particles in human respiratory tract. In this approach, the fluid phase is modeled with Navier-Stokes equations, and the particle phase is modeled using a large amount of individual particles obeying the Newtonian laws of motion. Some studies focused on investigation of specific type particles. Zhang et al. [8] analyzed the transport and deposition of cigarette smoke particles in CT-scan human airway model. They employed SST k ω − turbulence model for the airflow. Augusto et al. [9] studied the deposition of pharmaceutical aerosols in human airways using ideal human airway model under different breathing conditions. Tian et al. [10] investigated deposition patterns of inhaled wood dust in nasal cavity using CFD techniques. Zho et al. [11] investigated the deposition of fiber particles using empirical model which was developed based on the fiber deposition data from two human lung casts.
The human right and left lungs are asymmetrical. The right lung has 3 lobes: right upper lobe (RUL), right middle lobe (RML) and right lower lobe (RLL). The left lung has 2 lobes: Left upper lobe (LUL) and left lower lobe (LLL). The dimensions of the airways in same generation in the left and right lung are not the same. Hence airflow and particle transport and deposition study in human tracheobronchial airways is more meaningful if we take asymmetric model of the airways.
Many studies were carried out to see the ventilatory effects of dust on human respiratory system [12][13][14][15][16][17]. Although there are a number of studies on particle tracking and deposition in human respiratory tract, the transport and deposition of dust particles has not reported widely. The objective of this study was to investigate the trajectories and deposition of inspired dust particles in human tracheobronchial airways using CFD simulation. The characteristics and pattern of airflow go along with particle motion were also intended to be analysed.
Since sources of dust particles commonly inhaled by humans in roads and working environments are minerals of density 2 to 3 g/cm 3 , the density of dust particles approximated to be 2.65g/cm 3 [18][19][20]. The aerodynamic diameter of dust particles ranges from 0.5 to 50 µm [21]. Particulate matter of aerodynamics diameter less than 10 µm can pass through larynx and penetrate into the tracheobronchial airways [22][23][24]. In this study, dust particles of mass median aerodynamic diameter 10 µm was taken for the simulation.

Airway Geometric Model and Mesh Generation
Asymmetric bifurcation tracheobronchial airway model having 19 outlets was constructed using SOLIDWORKS 2013, 3D CAD software. The length and diameter of the airways were based on the CT scan data of adult male taken by Choi [25]. The branching angles of Horsfield et al. [26] were adopted for bifurcations. The mathematical description of the morphologically realistic bifurcation model of Hegedus et al. [27] was taken into consideration in the construction to make smooth transitions between parent and daughter tubes. We divided the airway model into seven zones to analyse the deposition of the dust particles. The constructed model and the regional zones are shown in Figure 1.
After importing the 3D airway model in ANSYS 19.0 Workbench, unstructured tetrahedral meshes with inflation layers were generated. Mesh independence study on the airflow solution was performed by comparing the average velocity of ten points on axis of plane cross-section at AA' (see Figure 1b) of the left main bronchus using grids consisting of 430016, 693695 and 1335828 cells. Steady simulation was implemented at each set of grids with inlet flow rate 15 L/min at the trachea. The increase of grid size from 693695 to 1335828 cells altered the average velocity with relative error 0.006. Considering the available memory and computational time, a mesh with 1335828 cells was adopted for the simulation.

Airflow
Air within the human respiratory tract may be considered to be a homogeneous, Newtonian and incompressible fluid [28]. The governing equations for airflow in respiratory tract are the continuity and Navier-Stokes equations for viscous incompressible Newtonian fluid. In vector notation these are: (2) Dust Particles in Human Tracheobronchial Airways Taking U. S. Standard Atmosphere Air Properties Data at sea level, the density and dynamic viscosity of air are taken to be ρ = 1.225 kg/m 3 and µ = 1.7984x10 -5 kg/(ms).
For light activity breathing condition (flow rate 30L/min) [29], the maximum airflow Reynolds number in the airway model was 2318 (at the trachea). The SST k ω − turbulence model was chosen to simulate the transitional and laminar flow regimes. The SST k ω − model is one of the Reynoldsaveraged Navier-Stokes (RANS) turbulence models which is more accurate and reliable for a wider class of flows as discussed by Elcner et al. [30].

Dust Particles Transport
For particle phase, we use Lagrangian discrete phase model (motion of individual particles is determined in space and time). If a particle has mass p m and tracks with velocity p v , then from Newton Second Law of Motion, we have where i F 's s are forces acting on the particle. The forces acting on the dust particles considered to be the drag force (F D ), gravitational force (F G ) and buoyancy force (F B ) and dust particles deposition in the airway model is mainly caused by inertial impaction and gravitational sedimentation. The drag force can be expressed as [8].
where D C is the drag coefficient, p A is the effective crosssectional area of the particle, and U is the flow velocity relative to particle velocity given by where p ρ is the density of the particle. From Eqs. (3-6), we obtain If p d is the diameter of the particle, we get The particle Reynolds number Re p is defined as From Eq. (9) and Eq. (10), we get where f D (v -v p ) is the drag force per unit mass of the particle with 2 Re 3 4 The drag coefficient C D is to be taken from Morsi

Boundary Conditions
For the airflow phase, the inlet velocity 0 1.8012 U = m/s was set at the trachea corresponding to the light activity breathing condition. The pressure outlet (zero relative pressure) was set at all outlet bronchi and no-slip boundary condition was imposed on walls of the airway model. For discrete phase model conditions, dust particles were released from the inlet surface with velocity equal to the air velocity. Escape was set at the inlet and at all outlets and trap at the wall of the airway model. It was assumed that a dust particle is deposited when it contacts with the airway wall. At the trachea, the Stokes number for the flow can be calculated as [29,32]  which shows that the particle motion followed the airflow pattern at inlet of the trachea. Here D = 0.0188m is the diameter of the trachea in the model.

Numerical Methods
Numerical simulation was performed with the commercial software ANSYS Fluent 19.0. Steady state solver model was used in the simulation. The second order discretization scheme was used for the pressure term and the second order upwind discretization scheme for convection terms. Green-Gauss cell based method was used to interpolate the diffusion gradients. The SIMPLE algorithm was applied for the pressure-velocity coupling. A residual less than 5 10 − was used as the convergence criteria for airflow simulation. Discrete phase model (DPM) with steady particle tracking was used to investigate the transport and deposition of the dust particles in the airway mode. Surface injection type with rosin-rammler diameter distribution was selected to inject the dust particles at the inlet with the same velocity as airflow using face normal direction. To ensure deposition independence, different numbers of particles (6900, 11500, 18400 and 27600) were taken at the injection. We found no significant difference in deposition when the number of particles injected is increased from 18400 to 27600. Based on this, injection with 18400 particles was selected in the study.

Validation
For validation of the computational methodology we used for airflow, velocity profiles of our simulation result on the right main bronchus cross-section at BB' (see Figure 1b) were compared with the results measured by De Rochefort et al. [33] using magnetic resonance (MR) gas velocimetry on a human lung model under same boundary conditions. The comparison is displayed in Figure 2. It was observed that the cross-sectional area on the bronchus of our airway model is much smaller than that of [33]. The simulated velocity profiles well agreed with the measured ones. The discrepancies may be because of differences in geometries of the airway models.  Figure 3 shows the airflow velocity streamlines of the simulation in the airway model. The streamlines pattern was relatively more complicated in right bronchus which may be resulted from geometric structure of the airway and the nearby bifurfications. The maximum velocity reached was 4.11 m/s. The velocity at end bronchi was generally higher than that at the main airways because of the narrow structure which is consistent with [34]. Figure 4 illustrates the spatial distribution of wall pressure in airway model. The wall pressure decreased down to the bronchi [35]. The pressure drop in the airway model was 13.08 Pa. The wall pressure on the right main and inferior lobar bronchi was relatively higher than that on left main and inferior lobar bronchi.  The wall shear stress distribution in the airway model is depicted in Figure 5. It is observed that the variation of wall shear stress distribution matched with those of the airflow velocity [34]. The maximum wall shear stress reached was 0.82 Pa.

Airflow Velocity, Wall Pressure and Wall Shear Stress
Velocity(m/s) Figure 5. Wall shear stress distribution in the tracheobronchial airway model. Figure 6 illustrates tracks of 600 dust particles from the total of 18400 tracks coloured by the particles residence time. The maximum residence time for the particles in the airway model was 0.322s. The location of 600 dust particles at different time after injection is illustrated in Figure 7. Particles escape at last from RUL bronchi outlets. The particles trajectories followed the airflow streamlines and the influence of secondary flow was high on the motion of the particles as reported by [36]. The deposition fraction of dust particles in different regional zones (see Figure 1b) of the tracheobronchial model is presented in table 1 and Figure 8. Out of 18400 injected dust particles 11007particles escaped from the airway model and 7393 particles were tapped in it. Higher deposition fraction (14.79%) was observed in the carina of the trachea, main bronchi and intermediate bronchus (Zone 1). The deposition fraction was higher in the right bronchial trees than that of the left bronchial trees which agrees with the result predicted by [37]. This result showed that there is a possibility that the right bronchial trees to be more affected by dust than the left bronchial trees. There was no particle deposited in the trachea in the simulation because the airflow and particle injection were set into its inlet in parallel with the vertical orientation of the trachea [38]. The total deposition fraction of dust particles in the tracheobronchial airway model was 40.18%. This is high relative to the small portion of the surface area that it covers in lung.

Conclusions
In this study, the airflow characteristics and dust particles transport and deposition were investigated in asymmetric human tracheobronchial airway model using CFD simulation. Because of geometric difference of the right and left lungs, there were differences in airflow characteristics, dust particles motion and deposition in the lungs. These facts were observed in the numerical simulation on the bronchial trees. The dust particles trajectories followed the airflow streamlines in the simulation. The deposition fraction of dust particles in human tracheobronchial airways was high particularly at the carina of trachea and the main bronchi. More dust particles were deposited in the right bronchial trees than in the left. This shows that the likelihood of the right lung and the carina to be affected by dust pollution is high. The deposition fraction of dust particles in human tracheobronchial airways was high. This study can provide awareness on deposition of dust particles passing beyond the larynx and enhance prevention of their entry into the respiratory system. It can also contribute a convenient way on the location of deposition of particles of a given type in human respiratory tract to be used in respiratory care efforts.

Conflicts of Interest
The authors declare that they have no conflicts of interest.