CFD–DEM simulation of particle transport and deposition in pulmonary airway

Powder Technology 228 (2012) 309–318

Contents lists available at SciVerse ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

CFD–DEM simulation of particle transport and deposition in pulmonary airway Xiaole Chen a, Wenqi Zhong a,⁎, Xianguang Zhou b, Baosheng Jin a, Baobin Sun b a b

Ministry of Education of Key Laboratory of Energy Thermal Conversion and Control, School of Energy and Environment, Southeast University, Nanjing, 210096, China Medical School of Southeast University, Nanjing 210009, PR China

a r t i c l e

i n f o

Article history: Received 19 January 2012 Received in revised form 16 May 2012 Accepted 20 May 2012 Available online 26 May 2012 Keywords: Airway CFD–DEM Deposition efficiency Particle motion Gas-solid flow

a b s t r a c t Including the effects of particle–particle collision and particle rotation as well as the volume of particle occupied in the fluid, the Computational Fluid Dynamics (CFD)–Discrete Element Method (DEM) approach could properly predict the behavior of gas–solid flow with strong coupling and the motion of non-spherical particle, thus show good potential in simulations of alveolar region and fibrous particle transport and deposition. CFD–DEM has been developed to investigate the particle transport and deposition characteristics in human airway. An idealized pulmonary airway model of generations 3 to 5 has been established using the same geometric parameters as previous experiment conducted by Kim and Fisher (C. S. Kim and D. M. Fisher. 1999). The predicted deposition efficiencies are in good agreement with experimental data. Thus, CFD–DEM could properly simulate the particle behaviors in pulmonary airway. Based on the simulations, particle motions are studied by analyzing the particle positions at different time intervals. The results show that the initial position of particle would notably affect its trajectory. The near wall particles distribute evenly in the daughter tubes when they cross the airway. The most centered particles at inlet travel through the model via the inner side tubes of generation 5. The trajectories of other particles would shift from the inside tubes to the lateral ones of last generation as the distance between particle initial position and tube center increases. © 2012 Published by Elsevier B.V.

1. Introduction Inhalable particles, especially the PM10 which has been proved to have strong association with respiratory disease and related mortality [1–3], arouse the interest of two phase flow researchers world wide. The transport and deposition characteristics of these particles are of fundamental importance. A thorough understanding could benefit the design of drug delivery device and elevate the efficacy of inhaled medicine. Many experiments have been carried out and provide valuable data. Kim et al. [4] investigated aerosol deposition in a single bifurcation airway. Kim and Fisher [5] continued their deposition research in a double bifurcation airway. Grgic et al. [6] studied the flow pattern and particle deposition in idealized mouth and throat models. Grgic et al. [7] examined the effect of unsteady flow rate increase on in vitro mouth–throat particle deposition. Golshahi [8] measured micron-particle deposition in nasal cavities of children and adolescents. Liu et al. [9] also studied particle deposition with their proposed “Carleton-Civic” standardized geometry of human nasal cavity. Chhabra and Prasad [10] reported airflow patterns and particle transport in an enlarged elastic alveoli model with particle image

⁎ Corresponding author at: Sipailou 2#, Nanjing 210096, Jiangsu, PR China. Tel.: + 86 25 83794744; fax: + 86 25 83795508. E-mail addresses: [email protected] (X. Chen), [email protected] (W. Zhong), [email protected] (X. Zhou), [email protected] (B. Jin), [email protected] (B. Sun). 0032-5910/$ – see front matter © 2012 Published by Elsevier B.V. doi:10.1016/j.powtec.2012.05.041

velocimetry (PIV) measurement. Bauer et al. [11] measured flow partition in a transparent silicone upper human lung airway model with oscillating flow. However, both the nasal cavity and pulmonary airway have complex geometry structures and significant individual differences. Thus, the fluid flow and particle transportation measurements as well as large sample experiments are always thorny subjects. The numerical simulation provides a possible solution to these problems. Detailed airflow and particle behavior could be obtained with proper gas–solid two phase flow simulation method. In the last two decades, many numerical investigations have been performed to study the aerosol deposition in human airway. Comer et al. [12] simulated the particle transport and deposition in airways with sequentially bifurcating, achieved good agreement with experimental results and analyzed the initial position of deposited particles. Then they [13] studied the air flow fields in the Weibel's G3–G5 airways, and discussed the differences of flow structures, energy dissipation rates and pressure drops in different models. They [14] also numerically researched the particle transport and deposition in these airways, and indicated that the particle release position significantly influence the deposition. Many other elements affecting particle deposition patterns and efficiencies in lung airways were analyzed in literature, i.e., the asymmetric branch flow rates [15], particle inlet distributions [16], curved inlet tubes [17], transient respirations [18–20], tumors [21] and nanoparticle [22,23]. Zhang and Kleinstreuer [24] compared the performance of four different turbulence models, i.e. the renormalization group (RNG) k-ε, Menter k-ω, Low-Reynolds-number (LRN) k-ε and LRN k-ω

310

X. Chen et al. / Powder Technology 228 (2012) 309–318

models, for simulating laminar–transitional–turbulent flows in constricted tubes and concluded that the LRN k-ω model was more suitable for such flow behavior simulations in constricted tubes than the other three methods. Kleinstreuer et al. [25] investigated gravitational effect on microparticle deposition in small airway models and indicated that sedimentation may change deposition pattern in deep airways especially under slow inhalation conditions. Kleinstreuer and Zhang [26] proposed a new methodology to compute particle deposition in large segments of human lung airways and compared deposition patterns and concentrations for nanoparticles and micron particles. Xi and Longest [27] evaluated the aerosol transport and deposition in nasal cavities based on a novel drift flux model with a near-wall velocity correction. Liu et al. [28] indicated that large eddy simulation (LES) method could provide more accurate particle deposition simulation results in nasal cavities than Reynolds averaged Navier–Stokes (RANS) with Lagrangian stochastic eddy interaction-models (EIMs) turbulent tracking and mean flow tracking. Moreover, the optimal drug-aerosol delivery, which could facilitate reducing side effect and dosage of expensive drug, drew increasing attention in recent years [29–32]. However, previous researches are mainly focused on the pulmonary airway variation and effect of breathing condition difference on the statistical results of particle deposition. The particle distribution and motion during the transport in human airway have not been clearly revealed. Targeting the rock mechanics at first, the general method of discrete element method (DEM) was invented by Cundall and Strack [33] and firstly known as the distinct element method. The DEM soon developed into several branches: the generalized discrete element method [34], the discontinuous deformation analysis [35] and the finite–discrete element method [36,37]. The DEM follows Newton's second law to calculate the individual element's motion. The DEM could consider friction, contact plasticity, gravity, cohesion and other interactions, record the information of each individual particle during simulation, e.g., position, velocity, forces exerted on, etc. Thus, it is widely applied in agriculture, pharmaceutical, mining, mineral processing and other researches. The rotation effect of spherical particle on the transport and deposition is limited and therefore is commonly not included in the simulations. However, the circumstance of non-spherical particle might not be the same. The experimental [38–42] and simulated [43–45] results both indicate that inertial impaction is the dominant mechanism of fibrous particle deposition. Thus, the rotation of non-spherical particle, which could lead to the change of particle's drag force during the transport, would be important in the non-spherical particle investigation in airways. The DEM provides a possible solution for the construction and motion description of irregular shaped particles using the multisphere method [46–50]. Considering the volume of the particle in fluid phase, the CFD–DEM could properly predict gas–solid flow with strong coupling [51–55]. The particle volume cannot be ignored in the pulmonary alveoli region simulations since the size of alveolus is very limited. Thus this method could be utilized to predict particle transport and deposition in these simulations. The possible applications above are also the two of many future focuses indicated by Kleinstreuer and Zhang [56]. However, the CFD–DEM has never been applied to biological multiphase flow before. Thus, it is important to examine whether these elements of CFD–DEM could properly simulate aerosol deposition in respiratory system firstly. The aim of this study is to evaluate whether the CFD–DEM is suitable for multiphase flow prediction of human respiration. Once this is validated, it would be interesting to investigate the particle distribution and motion in the airway during respiration. In order to compare the numerical method with experimental data, an idealized 3–5 generation pulmonary airway model has been constructed and shared the same geometry with previous experiment [5] which has been widely used

to validate simulated results [12,18,25]. The predicted deposition is compared with experimental data and the particle distributions and positions of different simulation time are also analyzed. 2. Methodology 2.1. Governing equations The fluid phase is considered as incompressible laminar flow. The continuity and the conservation of momentum equations are defined as
→ ∂ðερÞ=∂t þ ∇⋅ ερ u ¼ 0

ð1Þ


→
→→ 

→
→tr  → → þ ρε g − F ∂ ερ u =∂t þ ∇⋅ ερ u u ¼ −∇p þ ∇⋅ εμ ∇ u þ ∇ u

ð2Þ →

where ρ is the gas density, ε is the local void fraction, uis the vector of → gas velocity, p is the static pressure, g is the gravitational force,
→tr → → ∇ u is the transpose of ∇ u and F is the volumetric particle–fluid interaction force given by →

F¼

n → X fD;i =ΔV

ð3Þ

1 →

where fD;i is the drag force exerted on particle i, n is the number of particles in the specific computational cell and ΔV is the volume of the cell. The parabolic air inlet and other configurations of gas phase are explained in our previous CFD–DPM research [57]. The motion of particle is governed by Newton's second law → →

→

→

mp;i du p;i =dt ¼ mi g þ fD;i þ fc;i →

ð4Þ

→

Ii dωp;i =dt ¼ ∑Ti

ð5Þ →

→

where mp,i is the mass of particle i, u p;i and ωp;i are the translational →

and angular velocity vectors of particle i, fp;i is the contact force →

exerted on particle i, Ii is the particle i's rotational inertia and ∑Ti is the vector sums of torques caused by contacts of the particle. The model of DEM above is referred as model B in the literature [58–61]. The drag force is commonly the domain force driving the particles in airway research and is given as → →     u−u p;i =8


→ →  2 FD ¼ πdp C D u−u p;i

→

ð6Þ

where dp is the particle diameter, and CD, the drag coefficient is defined as 2

C D ¼ a1 =ReN þ a2 =ReN þ a3

ð7Þ

where  ReN is the Reynolds number for spherical particle → →  ReN ¼ ρ u −u p;i Dp =μ, a1, a2, and a3 are ReN determined coefficients [62]. Gravity is the natural field force exerted on particles. However, the gravitational force was sometimes neglected in previous upper airway simulations [63,64] because of the high velocity of inhaled air and short acting time of gravity. But Kleinstreuer et al. [25] indicated that gravity had significant influence on particle deposition in 8–9 generation lung airway. Since our pulmonary model has 3–5 generations, which was the middle section of the airway, the gravitational force was included in the calculation.

X. Chen et al. / Powder Technology 228 (2012) 309–318

f cn ¼

D3

qffiffiffiffiffiffiffi 3 4 Eeq Req δ2n 3

D3

Hertz–Mindlin and Deresiewicz contact theory [65] modified with damping force [66] was applied for modeling the particle–particle collision. The normal and tangential displacements are illustrated in Fig. 1 and the detailed equations are given as below. The normal elastic force fcn

311

ð8Þ

Generation 5 Generation 4

D2

where δn is the normal displacement. The equivalent Young's modulus Eeq and equivalent radius Req above are defined as

 2 2 1=Eeq ¼ 1−υi =Ei þ 1−υj =Ej

ð9Þ

1=Req ¼ 1=Ri þ 1=Rj

D1

ð10Þ

where Ei,υi, Ri and Ej,υj, Rj are Young's modulus, Poisson ratio and radius for element i and j, respectively. d The normal damping force fcn d f cn

Fig. 2. Schematic diagram of sequentially bifurcating model geometry.

d The tangential damping force fct

rffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5 rel ¼2 Sn meq U n lne= ln2 e þ π 2 6

ð11Þ d

where e is the coefficient of restitution, Unrel is the normal component of the relative velocity, meq is the equivalent mass and Sn is the normal stiffness defined as 1=meq ¼ 1=mi þ 1=mj Sn ¼ 2Eeq

ð12Þ

qffiffiffiffiffiffiffiffiffiffiffiffi Req δn

ð13Þ

with mi, mj the mass of particle elements i and j, respectively. The tangential elastic force fct  f ct ¼

St δt μ p f cn

f ct bμ p f cn f ct ≥μ p f cn

ð14Þ

with δt as the tangential displacement, μ p as the coefficient of static friction and St as the tangential stiffness which is given as St ¼ 8Geq

Generation 3

qffiffiffiffiffiffiffiffiffiffiffiffi Req δn

ð15Þ

where Geq is the equivalent shear modulus defined as
 1=Geq ¼ ð2−υi Þ=Gi þ 2−υj =Gj

ð16Þ

with Gi and Gj as the shear modulus of elements i and j, respectively.

δn

fct

f ct ¼ 2

rffiffiffiqffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5 rel St meq U t lne= ln2 e þ π2 6

ð17Þ

where Utrel is the tangential component of the relative velocity. The human airway is covered by liquid which consists of the viscous mucus layer and watery periciliary liquid (PCL). Foreign particles firstly deposit on the mucus layer and then they are transported and extruded with mucus layer by PCL under the cilia beating effect, namely the mucociliary clearance. Since the PCL does not affect the particle deposition and the mucus transport velocity, which is 2.4– 5.5 mm/min [67], comparatively much slower than the velocities of inhaled air and particles, the contact between particles and airway liquid has been simplified that particles would deposit once they collide with airway boundary in this numerical work and the boundary is assumed to be very soft. 2.2. Model geometry In order to evaluate the CFD–DEM simulating particle transport and deposition in human airway, a 3–5 generation pulmonary model was generated based on Weibel's 23-generation lung model [68] with 60° bifurcation angle. It shares the same geometric parameters with Kim and Fisher's experimental work [5]. The tubes of generations 3 and 5 had been prolonged to install the pipes in their experiment. The model geometry is shown in Fig. 2, in which D1, D2 and D3 stands for the deposition areas of prior tube, first and second bifurcations. The CAD package Gambit was used to construct the airway model and create the mesh of the domain. The cell number of the three generation model was 2 172 414. Detailed parameters are summarized in Table 1. 2.3. Calculation configurations In order to capture the particle–particle and particle–wall collisions, the calculations of solid and fluid phases were separated. The time step

δt fcn

Table 1 Geometric parameters of airway model. Parameters

Fig. 1. Sketch of normal and tangential displacements.

Diameter/mm Length/mm

Generation number 3 6 100

4 5 18

5 3.5 19

312

X. Chen et al. / Powder Technology 228 (2012) 309–318

Table 2 Physical and numerical parameters. Parameters

Value

dp, particle diameter/μm ρp, particle density/kg/m3 Gp, shear modulus of particle/Pa Gw, shear modulus of boundary/Pa υp, Possion's ratio of particle υp, Possion's ratio of boudnary e, coefficient of restitution μp, friction coefficient of particle to particle μp, friction coefficient of particle to wall TR, Rayleigh time/s

10 125, 250, 375, 500, 625 2 × 105 5 × 105 0.3 0.3 0.9 0.3 0.7 4.24×10−7, 6.00×10− 7, 7.35×10−7, 8.49×10−7, 9.49×10− 7 0.0001 2.12×10−7, 3.00×10− 7, 3.68×10−7, 4.25×10−7, 4.75×10− 7

time step for fluid phase/s time step for solid phase/s

for particles is determined by Rayleigh time and is greatly smaller than fluid phase. The Rayleigh time TR is the propagation time of the shear wave traveling through a solid particle which is given by Ning [69]

T R ¼ πRp

ρp Gp

!1 2


 = 0:1631υp þ 0:8766

ð18Þ

with Rp, ρp, Gp, and υp as the radius, density, shear modulus and Poisson's ratio of the particle, respectively. The time steps of solid and fluid phases were set to be 50% Rayleigh time and 0.0001 s. Before the particles were put into the computational domain, a parabolic air inlet with an average velocity of 4 m/s was set and the fluid phase was calculated firstly until convergence (dimensional residuals less than 10 − 4). Afterwards, 10,000 particles were generated randomly near the model inlet. The particles were carried by air stream and then deposited on the boundary or flowed out of the control volume. The particle diameters used in Kim and Fisher's experiment were 3–7 μm [5]. However, the DEM simulation requires intensive computational resources. Thus, we assumed that the closely spaced particles can be viewed as a cluster of 10 μm diameter spherical element in order to reduce run time. This method has been validated in the previous CFD–DEM research of our group [70] and could successfully capture gas–solid interaction characteristics. Particle density was chosen as the argument to determine the Stokes number (St= ρddp2U/18μf D, μ f is the absolute viscosity of air and D is the diameter of generation 3), i.e. particle densities 125, 250, 375, 500 and 625 kg/m 3 for Stokes numbers 0.025873, 0.051745, 0.077618, 0.10349 and 0.12936. The parameters used in the simulation are summarized in Table 2. FLUENT was used to carry out numerical solutions for the fluid phase equations with user-defined functions, and EDEM was employed to calculate the solid phase motions. After each calculation, the collision and particle position data were drawn and analyzed by programs written by current authors. Thus, the quantity of deposited particles in each bifurcation areas could be determined. The simulations were carried out on a HP Z800 workstation with 8 × 3.20 GHz cores and 48 GB RAM. The average computing time for fluid flow and particle transport simulation was approximately 7 days.

Fig. 3. Experimental and predicted DEs in sequentially bifurcating pulmonary airway model.1The form of DE = 100(1 − 1/(aStb + 1)) [5] is used for fitted curve equation.

where mpd and mpz are the mass of deposited and entered particle in specified control volume. Due to the possible fluctuation of the flow rate and aerosol diameter, the Reynolds numbers varied during the experiments. In addition, differences between simulation and experiment may result from subtle variance in geometric model and deposition measurements. However, the simulation results agree well with experimental data generally. In order to evaluate the accuracy of the CFD–DEM approach, the data of fitted curves of the CFD–DPM [57] and CFD–DEM predicted results are compared with Kim and Fisher's experiment [5] and presented in Table 3. The simulations of CFD–DPM and CFD–DEM share the same geometric parameters, transient flow and other configurations. The average absolute differences of simulated and experimental DEs indicate that the current accuracy of CFD–DEM is at the same level with CFD–DPM. This could result from the dilute particle flow and the particle–cluster diameter assumption. However, the CFD–DEM may become an advantageous method in particle transport and deposition research in airway considering its aforementioned potentials in fibrous particle and alveolus simulations. 3.2. Particle motion Fig. 4(a) and (b) shows the x–z and x–y plane particle positions at different time intervals. Under the entrainment of the parabolic velocity air inlet [57], the particles also obtain similar velocity distribution and thus central particles travel faster than the ones near the boundary (Fig. 4(a) 0.00100032 s). In Fig. 4(a) 0.01850640 s, the particles which crossed generation 4 all enter the inside tubes of generation 5 while the particles which came into the inside and outside tubes are close in quantity in Fig. 4(a) 0.02300780 s. When the slower particles

Table 3 The deposition efficiencies of fitted curve of DEM, DPM simulations and experiment. Stokes number

3. Results 3.1. Deposition efficiency Fig. 3 demonstrates the measured and CFD–DEM predicted deposition efficiencies (DEs) at D1, D2 and D3 bifurcation adjacent regions, respectively. The DE is defined as DE ¼ mpd =mpz

ð19Þ

0.025 0.044 0.063 0.082 0.101 Mean absolute difference

Deposition efficiencies % First bifurcation

Second bifurcation

DEM

DPM

Experiment

DEM

DPM

Experiment

2.422 5.105 8.080 11.196 14.360 2.969

0.274 1.387 3.811 7.815 13.395 0.985

1.522 4.142 7.665 11.82 16.379 –

1.377 3.320 5.733 8.468 11.418 0.326

0.374 1.495 3.554 6.602 10.580 1.216

0.422 1.649 3.858 7.0794 11.227 –

X. Chen et al. / Powder Technology 228 (2012) 309–318

313

Fig. 4. The x–z and x–y planes' particle positions of different time: (a) 0.00100032 s to 0.04500790 s; (b) 0.220008000 s to 0.28000800 s (particle size has been enlarged to provide clear view).

cross generation 4, more particles flow into the outside daughter tubes. Fig. 4(b) depicts the behaviors of the slow particles, the initial positions of which site near the geometric boundary of the inlet. These particles distribute more evenly when they go through the bronchia. It can be concluded that particles behave differently depending on their initial positions which can be divided into two types. The near wall particles would disperse into the daughter tubes as they move across the airway. The trajectories of the rest near center particles are determined by their initial distances from tube center. The most centered particles would cross the airway through the inside bronchia and the trajectories would shift to the lateral tubes as the distances between particle initial positions and tube center increase. In order to analyze the particle behaviors quantitatively, the airway model is separated into several zones shown in Fig. 5 and the statistical data of particles are drawn from DEM results of different simulation time. Four zones' statistics, the end of G3 and three front halves of daughter tubes G4 and G5 were summarized in Fig. 6. The behaviors of centered fast particles and slow particles close to the boundary are depicted in Fig. 6(a)–(d) and (e)–(f), respectively. The detailed three dimensional vector averages of distance from tube center to particle and data of other zones were provided in the Appendix section.

-6 -4

6

-5 -3 -9 -2

5 9 3 2

4

1

-1 8 Zone Partition

± i Zone Number Tube Center Line

0 z 7 x

Fig. 5. Zone partitions and tube center line in airway model.

314

X. Chen et al. / Powder Technology 228 (2012) 309–318

a

b

c

d

e

f

Fig. 6. Number of particles in zone, distance of vector average vR and average distance R from tube center to particle of different times: (a) zone 0 of 0.001–0.0406 s; (b) zone 1 of 0.001–0.0406 s; (c) zone 3 of 0.001–0.0406 s; (d) zone 5 of 0.001–0.0406 s; (e) zone 0 of 0.220–0.280 s; (f) zone 1 of 0.220–0.280sxi yi zi The axis of local coordinates.

Three parameters, i.e. number of particle in zone, distance of vector average vR and average distance R from tube center to particle, are presented in Fig. 6. The distance of vector average vR is defined as vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !2 !2 !2ffi u n n n u X→ X X → → t xj =n þ yj =n þ zj =n vRi ¼ 1

1

ð20Þ

1

where vRi represents the distance of vector average of zone i, n is the → → → number of particle in zone i and the xj , yj and zj are the three dimensional vectors of particle j from tube center to particle, respectively. This parameter could indicate whether the particles are uniformly distributed in the tube. Once the locations of the particles are biased toward one side of the tube, the vR would increase. The average distance R is defined as n qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X → →2 → Ri ¼ xj þ yj2 þ zj2 =n

In Fig. 6(a), the distance of vector average of zone 0 remains at a low level. This indicated that the particles were located evenly around the tube center. However, the average distance has a logarithmic

z-5 y-5

x-5

z3 y3 z1

y1

x3

y3 x3

x1

z3

y3 ð21Þ

z

x3

1

with Ri as the average distance from tube center to particle of zone i.

x Fig. A.1. Schema of local Cartesian coordinates calculating vectors in different zones.

X. Chen et al. / Powder Technology 228 (2012) 309–318

a

315

Time 0.0010 0.0054 1.75

Average Position

0.0098

1.75 0.75 1.25

0.75 1.25

0.0142 0.0186 0.0230

Zone 3 and 4

Zone 5 and 6

0.0274 0.0318 0.0362 0.0406 Zone 0, 2, 3, 6

Zone 7, 1, 3, 5

2.5 2.0

1.0

Viewing Direction

z x

1.0

2.0

Zone 1 and 2

0.2

3.0 0.05

0.10

0.15

Zone 7 and 0 Fig. A.2. The shift of average position vectors of different zones: (a) 0.001–0.0406 s; (b) 0.220–0.280 s.

increase with time which suggests that the central particles leave this zone and thus have higher velocities. The average distance from tube center to particle gradually reduces from 2.208 mm to 1.758 mm after 0.0142 s in Fig. 6(b) and the average position vectors move towards the outer side in Fig. Fig. A.2(a). Thus it could be concluded that the particles firstly entered zone 1 leaning to the inner side of generation 4 and diffused (Fig. 4(a) 0.01850640 s). The following particles progressively shift to the lateral side of the tube and therefore the distance of vector average increases with their entering the zone. The number of particle has a stepwise decrease which could result from the reduction of particle concentration. Fig. 6(c) and (d) shows the change of particle number, which could provide a stronger evidence of the shift of the particle trajectories. The particle number of zone 3, i.e., the lateral tube of generation 5, has swiftly increased with the rapid decrease of particle number of zone 5. The increases of zone 3's distance of vector average in Fig. 6(c) and the

change of average position vectors imply that particles move towards the outside of the tube with time, which also could support the trajectory shift assumption. The characteristics of slower particles in zone 0 and 1 are plotted in Fig. 6(e)–(f). Since the boundary particles at inlet are much less than the central ones, the numbers of particles in zones 3 and 5 are limited and therefore these data are only supplied in the Appendix section. The average distance of zone 0 remains at 2.9 mm in Fig. 6(e), which suggests that those particles are close to the boundary due to the fact that the radius of generation 3 is 3 mm. As the particles gradually flow out of this region, particle number reduces, and thus the asymmetry of particle position rises, leading to the increase of distance of vector average. The low level of zone 1's average distance and distance of vector average in Fig. 6(f) compared to Fig. 6(b) implies that the distribution of particles is comparatively more uniform and thus those particles would enter both inside and outside the tubes of generation 5.

316

X. Chen et al. / Powder Technology 228 (2012) 309–318

b

Time 0.2200 0.2260 1.75

Average Position

0.2320

1.75 0.75 1.25

0.75 1.25

0.2380 0.2440 0.2500

Zone 3 and 4

Zone 5 and 6

0.2560 0.2620 0.2680 0.2740 0.2800 Zone 0, 2, 3, 6

Zone 7, 1, 3, 5

2.5 1.0

2.0

Viewing Direction

z x

1.0

2.0

Zone 1 and 2

3.0

Zone 7 and 0 Fig. A2 (continued).

4. Conclusion The CFD–DEM approach has been applied to gas–solid two phase flow in airway for the first time. The previous researches have proven that the CFD–DEM is a useful tool to construct nonspherical particle and simulates the strong coupled gas–solid flow because it includes the particle–particle collision, particle rotation, the particle occupied volume in fluid and other phenomena. Thus, this method shows good potential in the prediction of fibrous particle transport and deposition and simulation of air-particle flow in alveolar region. In this paper, the results of numerical modeling of spherical particle transport and deposition in three-generation pulmonary airway are presented. The DEM technique has been coupled with CFD calculated fluid phase to track the particle trajectories and collisions. The following conclusions can be drawn based on current simulation.

(1) The CFD–DEM predicted deposition efficiencies agree well with experimental results, and therefore it is suitable for the prediction of multiphase flow in human airway. (2) The particle trajectories of the 3–5 generations airway based on Weibel's model with parabolic air inlet depend on their initial positions. A. The most centered particles flow into the inside tubes of generation 5. B. The near wall particles disperse into the daughter tubes evenly when crossing the airway. C. The trajectories of other particles in generation 5 would shift from the inside to the lateral tubes as the distance between particle initial positions and tube center increase. The DEM could provide more detailed information about particle motion and distribution and thus it can facilitate the understanding of the microcosmic mechanism of respiratory multiphase flow.

X. Chen et al. / Powder Technology 228 (2012) 309–318

317

Table A.1 The average distance from tube center in different zones. Time s

average distance mm

0.0010 0.0054 0.0098 0.0142 0.0186 0.0230 0.0274 0.0318 0.0362 0.0406 0.0450 0.2200 0.2260 0.2320 0.2380 0.2440 0.2500 0.2560 0.2620 0.2680 0.2740 0.2800

Zone 0

Zone 1

Zone 2

Zone 3

Zone 4

Zone 5

Zone 6

Zone 7

Zone 8

Zone 9

– – 0.313465 1.40406 1.91844 2.16331 2.33313 2.45578 2.48125 2.5903 2.62322 2.92081 2.92298 2.92296 2.92611 2.92886 2.93045 2.93072 2.93181 2.93149 2.93853 2.95039

– – – 2.20811 2.06799 2.08967 2.03887 2.05342 1.94221 1.76861 1.79482 1.07744 1.16541 1.22956 1.36403 1.38534 1.50993 1.49184 1.41462 1.36077 1.26714 1.27229

– – – 2.27667 2.20167 2.18164 2.21263 2.08771 2.19221 1.84688 1.80841 1.22138 1.33583 1.25075 1.12933 1.09348 1.00046 1.0738 0.873823 1.00857 1.07358 0.933878

– – – – – 1.49173 1.35461 1.27625 1.2782 1.3457 1.18278 – 0.760516 0.655963 0.723431 0.566375 0.55498 0.563028 0.749507 0.674577 0.827026 0.866419

– – – – – 1.49369 1.34074 1.2918 1.31074 1.28848 1.19324 – 0.822461 0.848314 0.809789 0.777039 0.766641 0.753913 0.773651 0.886818 0.855613 0.749436

– – – – 1.33962 1.45723 1.42515 1.42917 1.45554 1.21278 1.2653 1.57345 0.753298 0.780114 0.845743 0.82053 0.893422 0.933643 0.858614 1.12547 0.963755 0.603407

– – – – 1.27504 1.33516 1.39266 1.41041 1.38936 1.54426 1.46288 – 1.08378 0.980382 1.51257 1.4178 1.3 0.93527 1.27642 1.28221 1.31974 1.57979

1.99416 1.99389 2.00006 2.34988 2.52656 2.6214 2.68966 2.73185 2.77608 2.78412 2.81738 2.93429 2.93458 2.93815 2.93988 2.94446 2.94529 2.94505 2.94655 2.94643 2.94577 2.94466

– – – 1.18286 1.83937 2.10075 2.30185 2.41273 2.54478 2.60051 2.6201 3.1522 3.1187 3.14659 3.16638 3.18806 3.17995 3.15667 3.1405 3.1392 3.19397 2.85174

– – – – 2.15104 2.22134 2.00477 2.35184 1.93672 1.83386 1.36563 0.96549 0.727371 0.850249 0.871619 0.693826 0.935381 0.802679 0.754908 0.853415 0.716142 0.601776

Acknowledgments The authors gratefully acknowledge the financial support of the Foundation for the National Natural Science Foundation of China (No. 51176035), Author of National Excellent Doctoral Dissertation of PR China (No. 201040), and the Program for New Century Excellent Talents of the Ministry of Education of China (No. NCET090300). Appendix A In order to provide a direct view, the vectors in a specific zone are calculated by a local Cartesian coordinates using the tube center line as its z axis and the schema is illustrated in Fig. A.1 and the shift of

Table A.2 The particle number in different zones. Time s

0.0010 0.0054 0.0098 0.0142 0.0186 0.0230 0.0274 0.0318 0.0362 0.0406 0.0450 0.2200 0.2260 0.2320 0.2380 0.2440 0.2500 0.2560 0.2620 0.2680 0.2740 0.2800

Particle number Zone 0

Zone 1

Zone 2

Zone 3

Zone 4

Zone 5

Zone 6

Zone 7

Zone 8

Zone 9

0 0 37 918 643 542 487 338 420 103 506 516 453 354 228 152 100 59 38 21 12 10

0 0 0 491 372 311 319 333 329 382 248 56 64 79 92 120 102 89 84 74 55 39

0 0 0 233 490 415 404 330 222 303 208 16 30 26 42 30 48 43 25 36 32 31

0 0 0 0 0 187 214 183 192 109 115 0 9 11 6 12 6 8 5 5 7 11

0 0 0 0 0 16 273 323 266 227 200 0 6 7 12 9 16 15 11 5 7 11

0 0 0 0 301 241 120 96 73 44 69 1 2 9 7 12 7 8 12 2 3 3

0 0 0 0 407 447 319 181 144 101 107 0 2 4 6 10 17 16 15 17 12 6

9999 9999 9962 7276 5560 4540 3757 3253 2766 2654 2201 192 139 116 90 69 65 64 59 58 55 55

0 0 0 390 292 198 180 179 181 78 23 87 121 136 164 126 84 68 46 30 17 6

0 0 0 0 141 135 75 55 52 37 55 4 4 10 10 18 16 14 12 8 12 5

average position vectors for different zones is plotted in Fig. A.2. The average distance from tube center to particles and particle number in different zones are presented in Table A.1 and Table A.2, respectively.

Reference [1] D.E. Abbey, N. Nishino, W.F. McDonnell, R.J. Burchette, S.F. Knutsen, W. Lawrence Beeson, J.X. Yang, Long-Term inhalable particles and other air pollutants related to mortality in nonsmokers, American Journal of Respiratory and Critical Care 159 (2) (1999) 373. [2] B. Peterson, A. Saxon, Global increases in allergic respiratory disease: the possible role of diesel exhaust particles, Annals of Allergy, Asthma & Immunology 77 (4) (1996) 263–270. [3] J. Schwartz, Short term fluctuations in air pollution and hospital admissions of the elderly for respiratory disease, Thorax 50 (5) (1995) 531–538. [4] C.S. Kim, D.M. Fisher, D.J. Lutz, T.R. Gerrity, Particle deposition in bifurcating airway models with varying airway geometry, Journal of Aerosol Science 25 (3) (1994) 567–581. [5] C.S. Kim, D.M. Fisher, Deposition characteristics of aerosol particles in sequentially bifurcating airway models, Aerosol Science and Technology 31 (2–3) (1999) 198–220. [6] B. Grgic, W.H. Finlay, A.F. Heenan, Regional aerosol deposition and flow measurements in an idealized mouth and throat, Journal of Aerosol Science 35 (1) (2004) 21–32. [7] B. Grgic, A.R. Martin, W.H. Finlay, The effect of unsteady flow rate increase on in vitro mouth–throat deposition of inhaled boluses, Journal of Aerosol Science 37 (10) (2006) 1222–1233. [8] L. Golshahi, M.L. Noga, R.B. Thompson, W.H. Finlay, In vitro deposition measurement of inhaled micrometer-sized particles in extrathoracic airways of children and adolescents during nose breathing, Journal of Aerosol Science 42 (7) (2011) 474–488. [9] Y.A. Liu, E.A. Matida, M.R. Johnson, Experimental measurements and computational modeling of aerosol deposition in the carleton-civic standardized human nasal cavity, Journal of Aerosol Science 41 (6) (2010) 569–586. [10] S. Chhabra, A.K. Prasad, Flow and particle dispersion in a pulmonary alveolus—part I: velocity measurements and convective particle transport, Journal of Biomechanical Engineering T Asme 132 (5) (2010). [11] K. Bauer, H. Chaves, C. Brucker, Visualizing flow partitioning in a model of the upper human lung airways, J Biomech Eng-T Asme. 132 (3) (2010). [12] J.K. Comer, C. Kleinstreuer, S. Hyun, C.S. Kim, Aerosol transport and deposition in sequentially bifurcating airways, Journal of Biomechanical Engineering Transactions ASME 122 (2) (2000) 152–158. [13] J.K. Comer, C. Kleinstreuer, Z. Zhang, Flow structures and particle deposition patterns in double-bifurcation airway models. Part 1. Air flow fields, Journal of Fluid Mechanics 435 (2001) 25–54. [14] J.K. Comer, C. Kleinstreuer, C.S. Kim, Flow structures and particle deposition patterns in double-bifurcation airway models. Part 2. Aerosol transport and deposition, Journal of Fluid Mechanics 435 (2001) 55–80. [15] Z. Zhang, C. Kleinstreuer, C.S. Kim, Effects of asymmetric branch flow rates on aerosol deposition in bifurcating airways, Journal of Medical Engineering & Technology 24 (5) (2000) 192–202.

318

X. Chen et al. / Powder Technology 228 (2012) 309–318

[16] Z. Zhang, C. Kleinstreuer, Effect of particle inlet distributions on deposition in a triple bifurcation lung airway model, Journal of Aerosol Medicine Depos Clear Effect Lung 14 (1) (2001) 13–29. [17] Z. Zhang, C. Kleinstreuer, C.S. Kim, Effects of curved inlet tubes on air flow and particle deposition in bifurcating lung models, Journal of Biomechanics 34 (5) (2001) 659–669. [18] Z. Zhang, C. Kleinstreuer, C.S. Kim, Gas–solid two-phase flow in a triple bifurcation lung airway model, International Journal of Multiphase Flow 28 (6) (2002) 1021–1046. [19] Z. Zhang, C. Kleinstreuer, C.S. Kim, Cyclic micron-size particle inhalation and deposition in a triple bifurcation lung airway model, Journal of Aerosol Science 33 (2) (2002) 257–281. [20] Z. Li, C. Kleinstreuer, Z. Zhang, Particle deposition in the human tracheobronchial airways due to transient inspiratory flow patterns, Journal of Aerosol Science 38 (6) (2007) 625–644. [21] Z. Zhang, C. Kleinstreuer, C.S. Kim, A.J. Hickey, Aerosol transport and deposition in a triple bifurcation bronchial airway model with local tumors, Inhalation Toxicology 14 (11) (2002) 1111–1133. [22] H. Shi, C. Kleinstreuer, Z. Zhang, C.S. Kim, Nanoparticle transport and deposition in bifurcating tubes with different inlet conditions, Physics of Fluids 16 (7) (2004) 2199–2213. [23] Z. Zhang, C. Kleinstreuer, Airflow structures and nano-particle deposition in a human upper airway model, Journal of Computational Physics 198 (1) (2004) 178–210. [24] Z. Zhang, C. Kleinstreuer, Low-Reynolds-number turbulent flows in locally constricted conduits: a comparison study, AIAA Journal 41 (5) (2003) 831–840. [25] C. Kleinstreuer, Z. Zhang, C.S. Kim, Combined inertial and gravitational deposition of microparticles in small model airways of a human respiratory system, Journal of Aerosol Science 38 (10) (2007) 1047–1061. [26] C. Kleinstreuer, Z. Zhang, An adjustable triple-bifurcation unit model for air-particle flow simulations in human tracheobronchial airways, Journal of Biomechanical Engineering T ASME 131 (2) (2009). [27] J.X. Xi, P.W. Longest, Numerical predictions of submicrometer aerosol deposition in the nasal cavity using a novel drift flux approach, International Journal of Heat and Mass Transfer 51 (23–24) (2008) 5562–5577. [28] Y. Liu, E.A. Matida, J. Gu, M.R. Johnson, Numerical simulation of aerosol deposition in a 3-d human nasal cavity using RANS, RANS/EIM, and LES, Journal of Aerosol Science 38 (7) (2007) 683–700. [29] C. Kleinstreuer, Z. Zhang, Targeted drug aerosol deposition analysis for a fourgeneration lung airway model with hemispherical tumors, Journal of Biomechanical Engineering Transactions ASME 125 (2) (2003) 197–206. [30] C. Kleinstreuer, Z. Zhang, J.F. Donohue, Targeted drug-aerosol delivery in the human respiratory system, Annual Review of Biomedical Engineering 10 (2008) 195–220. [31] C. Kleinstreuer, Z. Zhang, Z. Li, W.L. Roberts, C. Rojas, A new methodology for targeting drug-aerosols in the human respiratory system, International Journal of Heat and Mass Transfer 51 (23-24) (2008) 5578–5589. [32] C. Kleinstreuer, Z. Zhang, Optimal drug-aerosol delivery to predetermined lung sites, Journal of Heat Transfer-Transactions of the ASME 133 (1) (2011) 011002. [33] P.A. Cundall, O. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1) (1979) 47–65. [34] J.R. Williams, G. Hocking, G.G.W. Mustoe, The theoretical basis of the discrete element method, Numerical Methods in Engineering Theory and Application (1985) 897–906. [35] G. Shi, Discontinuous Deformation Analysis: A New Numerical Model for the Statics and Dynamics of Block Systems, University of California, Berkeley, 1988. [36] A. Munjiza, D. Owen, N. Bicanic, A combined finite-discrete element method in transient dynamics of fracturing solids, Engineering Computations 12 (2) (1995) 145–174. [37] N.M.O. ES, J. Dolbow, T. Belytschko, A finite element method for crack growth without remeshing, International Journal for Numerical Methods in Engineering 46 (1999) 131–150. [38] W.C. Su, Y.S. Cheng, Deposition of fiber in the human nasal airway, Aerosol Science and Technology 39 (9) (2005) 888–901. [39] W.C. Su, Y.S. Cheng, Deposition of fiber in a human airway replica, Journal of Aerosol Science 37 (11) (2006) 1429–1441. [40] W.C. Su, Y.S. Cheng, Fiber deposition pattern in two human respiratory tract replicas, Inhalation Toxicology 18 (10) (2006) 749–760. [41] Y. Zhou, W.C. Su, Y.S. Cheng, Fiber deposition in the tracheobronchial region: experimental measurements, Inhalation Toxicology 19 (13) (2007) 1071–1078. [42] Y. Zhou, W.C. Su, Y.S. Cheng, Fiber deposition in the tracheobronchial region: deposition equations, Inhalation Toxicology 20 (13) (2008) 1191–1198. [43] K. Inthavong, J. Wen, Z. Tian, J. Tu, Numerical study of fibre deposition in a human nasal cavity, Journal of Aerosol Science 39 (3) (2008) 253–265. [44] Z. Wang, P.K. Hopke, G. Ahmadi, Y.S. Cheng, P.A. Baron, Fibrous particle deposition in human nasal passage: the influence of particle length, flow rate, and geometry of nasal airway, Journal of Aerosol Science 39 (12) (2008) 1040–1054.

[45] K.T. Shanley, G. Ahmadi, P.K. Hopke, Y.S. Cheng, Fibrous and spherical particle transport and deposition in the human nasal airway: a computational fluid dynamics model, ASME 2009 Fluids Engineering Division Summer Meeting, 2009 FEDSM2009-78204:2139-2148. [46] J. Favier, M. Abbaspour-Fard, M. Kremmer, Modeling nonspherical particles using multisphere discrete elements, Journal of Engineering Mechanics 127 (10) (2001) 971–977. [47] L. Vu-Quoc, X. Zhang, O. Walton, A 3-D discrete-element method for dry granular flows of ellipsoidal particles, Computer Methods in Applied Mechanics and Engineering 187 (3) (2000) 483–528. [48] M. Abbaspour-Fard, Theoretical validation of a multi-sphere, discrete element model suitable for biomaterials handling simulation, Biosystems Engineering 88 (2) (2004) 153–161. [49] H. Abou-Chakra, J. Baxter, U. Tuzun, Three-dimensional particle shape descriptors for computer simulation of non-spherical particulate assemblies, Advanced Powder Technology 15 (1) (2004) 63–77. [50] M. Kodam, J. Curtis, B. Hancock, C. Wassgren, Discrete element method modeling of bi-convex pharmaceutical tablets: contact detection algorithms and validation, Chemical Engineering Science 69 (1) (2012) 587–601. [51] M.H. Zhang, K.W. Chu, F. Wei, A.B. Yu, A CFD–DEM study of the cluster behavior in riser and downer reactors, Powder Technology 184 (2) (2008) 151–165. [52] B. Ren, W. Zhong, B. Jin, Z. Yuan, Y. Lu, Modeling of gas‐particle turbulent flow in spout‐fluid bed by computational fluid dynamics with discrete element method, Chemical Engineering and Technology 34 (12) (2011) 2059–2068. [53] B. Ren, Y. Shao, W. Zhong, B. Jin, Z. Yuan, Y. Lu, Investigation of mixing behaviors in a spouted bed with different density particles using discrete element method, Powder Technology 222 (2012) 85–94. [54] H. Zhu, Z. Zhou, R. Yang, A. Yu, Discrete particle simulation of particulate systems: theoretical developments, Chemical Engineering Science 62 (13) (2007) 3378–3396. [55] C. Jayasundara, R. Yang, B. Guo, A. Yu, I. Govender, A. Mainza, A. Westhuizen, J. Rubenstein, CFD–DEM modelling of particle flow in IsaMills—comparison between simulations and PEPT measurements, Minerals Engineering 24 (3–4) (2011) 181–187. [56] C. Kleinstreuer, Z. Zhang, Airflow and particle transport in the human respiratory system, Annual Review of Fluid Mechanics 42 (2010) 301–334. [57] X. Chen, W. Zhong, B. Sun, B. Jin, X. Zhou, Study on gas/solid flow in an obstructed pulmonary airway with transient flow based on CFD–DPM approach, Powder Technology 217 (2012) 252–260. [58] D. Gidaspow, Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions, Academic Press, 1994. [59] K. Kafui, C. Thornton, M. Adams, Discrete particle-continuum fluid modelling of gas–solid fluidised beds, Chemical Engineering Science 57 (13) (2002) 2395–2410. [60] Y.Q. Feng, A.B. Yu, Assessment of model formulations in the discrete particle simulation of gas−solid flow, Industrial and Engineering Chemistry Research 43 (26) (2004) 8378–8390. [61] K.D. Kafui, C. Thornton, M.J. Adams, Reply to comments by Feng and Yu on “Discrete particle-continuum fluid modelling of gas–solid fluidised beds”, Chemical Engineering Science 59 (3) (2004) 723–725. [62] S.A. Morsi, A.J. Alexande, Investigation of particle trajectories in 2-phase flow systems, Journal of Fluid Mechanics 55 (26) (1972) 193–208. [63] Z. Zhang, C. Kleinstreuer, C.S. Kim, Micro-particle transport and deposition in a human oral airway model, Journal of Aerosol Science 33 (12) (2002) 1635–1652. [64] Z. Zhang, C. Kleinstreuer, J.F. Donohue, C.S. Kim, Comparison of micro- and nano-size particle depositions in a human upper airway model, Journal of Aerosol Science 36 (2) (2005) 211–233. [65] R.D. Mindlin, H. Deresiewicz, Elastic spheres in contact under varying oblique forces, Journal of Applied Mechanics-Transaction of the ASME 20 (3) (1953) 327–344. [66] A.O. Raji, Discrete element modelling of the deformation of bulk agricultural particulates, Newcastle upon Tyne, Newcastle University, 1999. [67] W. Foster, E. Langenback, E. Bergofsky, Lung mucociliary function in man: interdependence of bronchial and tracheal mucus transport velocities with lung clearance in bronchial asthma and healthy subjects, Annals of Occupational Hygiene 26 (2) (1982) 227. [68] E.R. Weibel, Morphometry of the Human Lung, Academic Press, New York, 1963. [69] Z. Ning. Elasto-plastic impact of fine particles and fragmentation of small agglomerates. PhD Thesis, Aston University, Birmingham. 1999. [70] B. Ren, W.Q. Zhong, B.S. Jin, Z.L. Yuan, Y. Lu, Computational fluid dynamics (CFD)–discrete element method (DEM) simulation of gas–solid turbulent flow in a cylindrical spouted bed with a conical base, Energy and Fuels 25 (9) (2011) 4095–4105.