A fully Eulerian approach to particle inertial deposition in a physiologically realistic bifurcation

Applied Mathematical Modelling 37 (2013) 5591–5605

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

A fully Eulerian approach to particle inertial deposition in a physiologically realistic bifurcation M. Pilou a,b,⇑, V. Antonopoulos b, E. Makris a,b, P. Neofytou a, S. Tsangaris b, C. Housiadas a a

Thermal Hydraulics & Multiphase Flow Laboratory, Institute of Nuclear Technology & Radiation Protection, National Centre for Scientific Research ‘‘Demokritos’’, 15310 Agia Paraskevi, Greece b Laboratory of Biofluid Mechanics & Biomedical Engineering, School of Mechanical Engineering, National Technical University of Athens, 15780 Zografou, Greece

a r t i c l e

i n f o

Article history: Received 23 April 2012 Received in revised form 18 September 2012 Accepted 15 October 2012 Available online 5 November 2012 Keywords: Particle deposition Physiologically realistic bifurcation Eulerian–Eulerian method Inertial effects CFD modelling 3D transport

a b s t r a c t The investigation of aerosol flows in bifurcations is significant in biomedical applications, because this geometry resembles the branching airways of the lung. Most of the existing computational studies regarding aerosol flow in branching airways use Lagrangian description for the particulate phase, when particle inertial effects are taken into consideration. In the present study a fully Eulerian model is used for studying the transport and deposition of inertial particles under steady state inspiratory flow in the single physiologically realistic bifurcation created by generations G3-G4 of the human lung. In-house codes are used to create the geometry and the computational grid, obtain the fluid flow field and calculate particle concentration. Deposition fractions are calculated, concentration profiles are shown and deposition sites are indicated under different flow conditions (Reynolds number and flow asymmetry in the branches) and particles of various sizes (1–10 lm). It is found that total particle deposition fraction is loosely dependent on fluid flow Reynolds number and our results are in agreement with experimental findings. Moreover, total deposition fraction does not change significantly between symmetric Q 1 =Q 2 ¼ 1 and asymmetric Q 1 =Q 2 ¼ 2 flow conditions, but is considerably lower for the totally obstructed case, Q 2 ¼ 0. On the contrary, deposition sites and particle concentration profiles depend on the various flow conditions. Apart from deposition at the bifurcation, which is present at all cases, deposition sites downstream the bifurcation move towards the outer bifurcation wall and daughter tube exit as Reynolds number increases. For the Q 2 ¼ 0, it is also shown that particles are trapped and deposit even in the obstructed daughter tube. In all studied cases, it is found that there is a particle-free region at the outer wall, at the beginning of the daughter tubes opposite the bifurcation. The straightforward formulation of deposition modelling and derivation of detailed concentration profiles along the bifurcation are major advantages of the fully Eulerian methodology used in the study. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction Aerosol flows in bifurcations are of great importance for biomedical applications, because this geometry can be considered as a building block of the lungs. In the literature, one can find experimental, analytical and computational studies regarding the transport and deposition of particles in single, double or multiple bifurcations. However, aerosol flows in these

⇑ Corresponding author at: Thermal Hydraulics & Multiphase Flow Laboratory, Institute of Nuclear Technology & Radiation Protection, National Centre for Scientific Research ‘‘Demokritos’’, 15310 Agia Paraskevi, Greece. Tel.: +30 210 6503700; fax: +30 210 6545496. E-mail address: [email protected] (M. Pilou). 0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.10.055

5592

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

Nomenclature Latin c Cc ~ cs dp d1 Fr ~ g Jc Jd ^ ^ı; ^|; k L ^ n p Pe Q Re ~ S St t ~ t u; m; w ~ vp

particle concentration Cunningham slip correction factor gravitational settling velocity vector particle diameter parent tube diameter Froude number gravity vector convective flux diffusive flux Cartesian basis vectors tube length normal to the surface unit vector fluid pressure Peclet number fluid flow rate Reynolds number surface area vector Stokes number time fluid velocity vector Cartesian components of fluid velocity particle velocity vector

Greeks

g lf qf qp sij sp dX

deposition fraction fluid dynamic viscosity fluid density particle density components of the stress tensor particle relaxation time elemental volume

Subscripts f fluid m mean p particle w wall Abbreviations PRB physiologically realistic bifurcation CDS central differences discretization CFD computational fluid dynamics

complex geometries depend on many different parameters and phenomena, thus generalized conclusions are not always reached [1]. There are experimental studies, where a hollow cast of the tracheobronchial tree in addition to one or more parts of the upper respiratory tract, such as the larynx and the nasal or the oral cavity, were used for measurement of particle deposition in the human respiratory cast. For example, Schlesinger et al. [2] measured deposition of particles with mean aerodynamic diameters between 2.5 and 8.1 lm in a hollow cast of the larynx and the tracheobronchial tree under steady inspiratory flow, and showed that these particles deposit preferably at bifurcations. Moreover they concluded that deposition presented a maximum in the third lung generation (G3). Later, a cast of the oral airway, which included the oral cavity, the pharynx, the larynx, the trachea and three bronchi generations, was used in order to study the deposition of particles (0.93–30 lm) under different flow rates [3]. It was found that impaction is the dominant mechanism, as deposition efficiency depends primarily on the Stokes number. However, most of the studies, either experimental or numerical, focus on a specific part of the respiratory tract, usually on two or three consecutive generations of the lungs, therefore in particle transport and deposition in a single or double bifurcation.

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

5593

Kim and Iglesias [4] studied experimentally deposition in a Y-shaped bifurcation of circular cross-section of monodispersed oil droplets of 3, 5 and 7 lm, varying the branching angle and symmetry as well as fluid flow Reynolds number (Re ¼ 1132—3397) and flow patterns in each of the branches. They found that deposition occurs near the bifurcation and increased with increasing Stokes number, whereas it is not affected significantly by branching angles lower than 45°, branching asymmetry and flow patterns in the daughter branches for the tested particles and Reynolds numbers. Moreover, they showed that bend models [5–7] cannot realistically describe the deposition in bifurcating airways, as might be assumed on the basis of the geometrical similarity of a bifurcation to a system of two adjacent bends. The aforementioned result triggered further experimental works on the aerosol flow in bifurcations. Deposition of monodispersed particle populations was measured in a symmetrical single bifurcation of circular cross-section under steadystate flow conditions [8]. Particle diameter, dp , air flow Reynolds number, Re, bifurcation angle and parent to daughter diameter ratio served as problem parameters. The study showed that inertial impaction was the most significant deposition mechanism for the particle sizes under investigation (dp P 1 lm) and the majority of particles deposit in a narrow region at the top of the carinal ridge. Moreover, Kim and coworkers [8] found that bifurcation angle, parent to daughter diameter ratio and flow Reynolds number may influence the local particle deposition patterns, but did not affect significantly the total deposition fraction. The latter depends only on the Stokes number and increases monotonically with it. Experimental studies in double bifurcations reached similar conclusions [9–11]. In addition, these studies concluded that deposition in the second bifurcation is generally lower than the one in the first, a fact that is attributed to the fluid flow downstream of the first bifurcation. Apart from the previously mentioned experimental works, there are also many computational studies of the fluid-particle flow in branching airways. All of them treat the fluid flow using an Eulerian formulation, but the majority of the studies adopt a Lagrangian approach for the particulate phase. Lee and Goo [12] simulated the inertial deposition of particles in a bifurcating channel of square cross-section, assuming equal parent and daughter sizes, and studied the effect of different characteristics of the geometry and the fluid flow, as well as the size of the particles. They showed that both the shape of the carinal ridge, rounded or sharp, and the flow Reynolds number in combination with the bifurcating angle influence the fluid flow field and the particle deposition. Moreover, in a similar study [13], the deposition in a single bifurcation of square cross-section based on the geometric characteristics of the third and forth generations of the lung (G3-G4) was studied. They took into account impaction and sedimentation of particles under steady flow for Reynolds numbers equal to 100 and 1000 and uniform and parabolic inlet fluid velocity profiles. Their results indicated that deposition due to impaction is sensitive in both fluid flow Reynolds number and inlet velocity profile, whereas deposition due to sedimentation is almost independent for particles bigger than 10 lm. More recently, simulation of aerosol flow in a single bifurcation of circular cross section, based again on the G3-G4 generations characteristics and for particles in the size range of 1–500 nm, showed that molecular diffusion is the dominant deposition mechanism for particles smaller than 10 nm, whereas interception cannot be neglected for high flow rates and even becomes the major mechanism for particles larger than 20 nm [14]. In the last decade, partly owing to great advances in computational power, the interest is shifted to the simulation of three or more generations of the tracheobronchial tree. Comer and coworkers [15] studied the three dimensional, steady, laminar, inspiratory aerosol flow in a symmetric physiologically realistic model of the third to fifth lung generation (G3-G5), i.e., a double bifurcation. Flow Reynolds number ranged between 500 and 2000 and particle Stokes number between 0.02 and 0.12. In addition, the effect of carinal ridge shape and planar or non-planar configuration of the second bifurcation were investigated. Their results showed that total deposition efficiency depends almost entirely on the Stokes number. Nevertheless, the shape of the bifurcation, inlet velocity profile and flow Reynolds number have an effect on the particle deposition patterns. Moreover, they showed that the particle concentration patterns are related to the corresponding secondary flows and that in the first bifurcation, most particles deposit because of inertial impaction. In contrast to earlier results [12], Comer and coworkers found that the shape of the carinal ridge may alter the deposition patterns to some extent, but affects little the total deposition efficiency in the bifurcations. Other numerical studies of double bifurcations are the ones by Longest and coworkers, who studied the aerosol flows of submicron particles (dp 6 1 lm) in models of generations G3-G5 and G7-G9, took into account turbulence and also simulated steady expiratory flow [16–18]. In addition, Zhang and coworkers studied fluid-particle flows in triple bifurcations (G0-G3 [19] and G3-G6 [20]). Finally, particle deposition have been studied in a five generations tree [21] and in a 16generations tree by combining adjustable triple bifurcation units [22]. All the aforementioned studies use a Lagrangian description of the particulate phase. An Eulerian formulation is adopted by a few researchers but only for the transport and deposition of submicron particles, for which inertial effects are not important [23–26]. A notable exception is the work of Gradon and Orlicki [27], who studied the transport of submicron particles in single PRBs created by generations 9–20 according to Weibel’s model of the lung [28]. In these generations inertia and diffusion are equally important particle deposition mechanisms. In particular, they solved simultaneously the equations of mass and momentum of both the air and the particles. In this work an Eulerian-Eulerian model [29] is used for studying the transport and deposition of inertial particles in a single physiologically realistic bifurcation (PRB) under steady state inspiratory flow. One-way coupling of the particulate phase is considered whereby the particle motion is affected by the air flow, but not vice versa. The effect of fluid flow (Reynolds number, asymmetry in the branches) and particle size is investigated. The fully Eulerian methodology used offers clear advantages:

5594

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

(a) it does not require time and effort consuming particle tracking histories in order to obtain accurate statistics; (b) it takes into account particle transport due to diffusion and inertia simultaneously, thus it is valid for a wide range of particle sizes; (c) it is fairly simple, as it decouples the mass and momentum equations of the particulate phase and the particle velocity is given only in terms of the fluid velocity and the its spatial derivative; and (d) particle concentration is directly calculated as natural part of the solution of the particle population balance equation (PBE). 2. Model description 2.1. Air flow ^ dS If ~ t ¼ u^ı þ m^| þ wk^ is the fluid velocity, p the pressure, qf the fluid density, sij the components of the stress tensor, d~ S¼n ^ the Cartesian basis ^ the normal to the surface unit vector, dX the elemental volume, and ^ı, ^|, k the elemental surface with n vectors, then the mass continuity equation in integral form is:

Z

~ t  d~S ¼ 0;

ð1Þ

S

and the momentum equations (Navier–Stokes) are:

Z Z Z Z d ^  d~ S; qf u dX þ qf u~ t  d~S ¼  p^ı  d~S þ ðsxx^ı þ syx^| þ szx kÞ dt X S S S Z Z Z Z d ^  d~ S; q m dX þ qf m~ t  d~S ¼  p^|  d~S þ ðsxy^ı þ syy^| þ szy kÞ dt X f S S S Z Z Z Z d ^  d~ S: q w dX þ qf w~ t  d~S ¼  p k^  d~S þ ðsxz^ı þ syz^| þ szz kÞ dt X f S S S

ð2Þ

The in-house CFD code developed by Neofytou and Tsangaris [30] is used to solve the Navier–Stokes equations for the threedimensional incompressible flow field of the carrier fluid. Incompressibility of air, which is the fluid of interest, is justified due to the low flow velocities in addition to the temperatures (about 37 °C) in the respiratory system. In particular, for the conditions in the cases of interest, the speed of sound is about 300 m/s, thus the Mach numbers are very low. The code uses a finite-volume methodology with a collocated arrangement of variables, while it enables multi-block computations. A pressure-correction equation is used and the coupling of velocity and pressure is dealt with using the SIMPLE algorithm [31]. In the code, the convection terms are discretized using the QUICK difference scheme, which is of third-order, whereas the central difference scheme (CDS) is used for the diffusion (viscous) terms. Herein, a solution is reached assuming steady-state, laminar flow of a Newtonian fluid of constant properties. The influence of the particulate phase on the air flow field is considered negligible (one-way coupling). A parabolic velocity profile at the inlet is chosen, and a no-slip fluid velocity condition is imposed at the wall boundaries. 2.2. Aerosol particles If c is the particle concentration and ~ v p the particle velocity, then the general dynamic equation for a monodispersed particle population in a flowing fluid is written as:

  @c þ r  c~ v p ¼ 0; @t

ð3Þ

where the internal processses, i.e., growth and coagulation/aggregation have been neglected. Eq. (3) is the mass conservation equation of the particulate phase, which under steady-state conditions becomes:





r  c~ v p ¼ 0:

ð4Þ

To incorporate particle inertial effects on Brownian diffusive transport, Pilou et al. [29] used a first order correction of particle velocity in particle relaxation time [32]. Taking also into consideration gravity, particle velocity can be approximated by

  ~ v p ¼ ~t  sp~t  r~t  Dr ln c þ O s2p ; 2

ð5Þ

where sp ¼ qp dp C c =18lf v and D ¼ kB T f =3plf dp the particle relaxation time and Stokes–Einstein diffusion coefficient. In these relationships, kB the Boltzmann’s constant, T f and lf the fluid temperature and dynamic viscosity, qp and dp the particle density and diameter, C c the Cunningham slip correction factor, and v the particle dynamic shape factor. The last two factors are unity for spherical particles of diameters greater than 1 lm, such as those considered in this study. The inertial term in Eq. (5) implies that the particle velocity field is compressible for inertial particles regardless the presence of diffusion and the incompressibility of the carrier fluid. The velocity due to gravitational settling, ~ cs , is equal to:

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

5595

!

~ cs ¼

qp d2p qf ~ 1 g; 18lf qp

ð6Þ

with ~ g the gravity vector [33], and for the aerosol under study, where qf  qp , becomes:

~ cs ¼

qp d2p ~ g ¼ sp~ g: 18lf

ð7Þ

Particle velocity is, thus, decomposed into a diffusive part, dependent only on the particle concentration gradient, and a convective part, ~ v c , that equals:

~ v c ¼ ~t þ sp ð~g  ~t  r~tÞ;

ð8Þ

and depends only on the fluid field and particle relaxation time. Therefore convective velocity (Eq. (8)) incorporates the inertial and gravitational effects of the particles. Introducing the approximation of the particle velocity into Eq. (4), we obtain:

r  ðc~ v c Þ ¼ r  ðDr cÞ:

ð9Þ

This equation is a modified steady-state convective diffusion equation, which in dimensionless and integral form is written as (to simplify notation we retain the same variables to denote dimensionless quantities):

Z

~ ¼ Pe1 c~ v c  dS

Z

S

~ rc  dS;

ð10Þ

S

where the dimensionless convective velocity is equal to:

  ~ v c ¼ ~t þ St Fr1~g  ~t  r~t :

ð11Þ

In the previous equations, St ¼ sp to =d1 is the particle Stokes number, with to the mean fluid velocity at the PRB inlet, Fr ¼ t2o =gd1 is the Froude number and Pe ¼ d1 to =D is the mass Peclet number, a function of the particle diameter through D. The parent tube diameter d1 is chosen to scale lengths, whereas to is used to scale the velocities. Eqs. 10,11, describe particle convection, inertia, gravitational settling and diffusion in an Eulerian formulation. Contrary to previous studies [27], the methodology adopted in the present work decouples the mass and momentum equations of the particulate phase, being thus simpler. The modified convection–diffusion equation,Eq. (10), is solved assuming (dimensionless) concentration of unity at the inlet of the PRB. Moreover, the totally absorbing wall boundary condition is imposed on the PRB walls:

cjw ¼ 0:

ð12Þ

However, the particle convective velocity just before the wall has a finite value because of the external forces, which gives a total deposition flux equal to:

   F dep ¼ F c w þ F d  ;

ð13Þ

w

 where F c w

8 > <0   F c w ¼ R > ~  : S c~ v c  dS

w

 and F  the diffusive flux:

 ~  6 0; if ~ v c  dS w  ~ ~ if v c  dS > 0;

ð14Þ

w

d

w

Z    ~  : F d  ¼ Pe1 rc  dS w

S

w

ð15Þ

The particle convective velocity in Eq. (14) is calculated at the computational grid point closest to the wall. Eq. (10) with the appropriate boundary conditions is solved numerically using an in-house CFD-based code [29]. The code calculates the convective velocity of the particles by the fluid velocity, through Eq. (11), and then proceeds to the solution of the modified convection–diffusion equation in three dimensions using a finite volume methodology. In the code, the convective term is discretized using a second-order deferred correction [34], while a second-order central difference scheme is preferred for the diffusive term. 2.3. Geometry and computational grid A physiologically realistic bifurcation (PRB) based on the geometric characteristics of the bifurcation created by the third and fourth lung generations (G3-G4) according to Weibel’s model [28] is used in the simulations (Table 1). These generations

5596

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605 Table 1 Geometrical characteristics of G3-G4 bifurcation according to Weibel [28]. G3-G4 bifurcation Parent diameter, d1

5:6  103 m

Parent length, L1

11  103 m

Daughter diameter, d2

4:5  103 m

Daughter length, L2

9:2  103 m 35°

Bifurcation angle, aprb

are chosen because they have been repeatedly used in literature [13,8,15], and it was experimentally shown that particles deposition peaks there [2]. The PRB geometry is constructed using the Design Modeler application of the commercial CFD package ANSYS 11 and user defined functions. In particular, the equations that describe a curved carinal ridge by Heistracher and Hofmann [35] are incorporated for the description of the bifurcation region. The centerlines of the parent and daughters tubes and the whole constructed surface of the PRB with the curved carinal ridge are shown in Fig. 1(a) and (b), respectively. The three dimensional multi-block structured grid of the PRB is generated using the methodology and in-house software developed by Makris and coworkers [36,37], originally developed for the generation of structured grids of patient-specific geometries obtained from medical images, such as abdominal aortic aneurysms or an abdominal aorta bifurcation. The advantage of this method is the preservation of the naturally complex biological geometries with a structured grid, which has lower numerical diffusion and converges faster compared to an unstructured one [38,36,39]. The basic steps of the method are:

Fig. 1. (a) Centerlines and geometrical characteristics of the G3-G4 bifurcation (due to geometrical symmetry only one daughter tube is shown), (b) constructed surface of the G3-G4 bifurcation with curved carinal ridge, (c) unstructured surface grid, and (d) structured volume grid and grid crosssection.

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

5597

(a) creation of the geometrical surface (analytically or by image reconstruction techniques using MIMICS or 3DSlicer); (b) generation of an unstructured grid over the surface; (c) division of the unstructured grid into open sub-grids (Pyformex or ANSYS ICEM CFD); (d) transformation of each of the unstructured sub-grids into structured grid; (e) reassembly of the structured sub-grids to form the initial surface and definition of grid blocks; and (f) generation of the three dimensional volume grid of the space included by the surface. In Fig. 1(c), the unstructured grid generated on the PRB surface is shown, whereas in Fig. 1(d) the structured volume grid as well as the grid cross-section are shown. A grid independence study indicated that a grid resolution of 9  105 nodes is sufficient in order to obtain the fluid flow and the particle concentration fields. 3. Results and discussion In the physiologically realistic bifurcation study, the developed fully Eulerian computational model is used to investigate the effect of fluid flow and particle size not only on the deposition fraction, but also on the concentration profiles and the particle deposition sites. The fluid and particle properties used in the simulations are shown in Table 2. Flow symmetry is assumed for the three different flow conditions shown, which correspond to inspiration of an adult man during rest (Re ¼ 464) and under light (Re ¼ 1132) and heavy (Re ¼ 1788) exercise. Moreover, for Re ¼ 1132 (light exercise) flow asymmetry is examined. If Q 1 ; Q 2 are the flows exiting daughter tubes 1 and 2, respectively (Fig. 2), then the cases of the flow in tube 1 being double than the one in tube 2, i.e., Q 1 =Q 2 ¼ 2, as well as tube 2 being totally obstructed, i.e., Q 2 ¼ 0, are investigated. The former corresponds to partial obstruction of the flow in the daughter tube 2, whereas the latter to totally obstructed flow. Fully developed fluid velocity (parabolic) profile and uniform (plug) particle concentration profile at the entrance of the third lung generation (inlet) are assumed for all cases. In Fig. 2, different cross-sections are defined. More specifically, cross-section A-A lies at the parent tube in the beginning of the bifurcation region (y = 2.4), whereas cross-sections B-B0 and C-C0 are located downstream the bifurcation, at y = 3.5 and y = 4.5, respectively. The subscripts 1 and 2 refer to the corresponding daughter tube. Note, however, that in the presentation of the results, the subscript is omitted in the symmetric flow cases. In addition, in the insert of Fig. 2, the diameters H-H and V-V of a cross-section are defined, the first of which lies parallel and the second perpendicular to the plane of the PRB geometric symmetry z ¼ 0. 3.1. Air flow The effect of Reynolds number to the axial velocity profiles along the diameters H-H parallel (top) to and V-V perpendicular (bottom) to the symmetry plane (z ¼ 0) for the symmetric flow case are shown in Fig. 3 at the different cross-sections. Upstream of the bifurcation (A-A cross-section, left), the axial velocity has a symmetric almost parabolic profile that is nearly independent of the Re along H-H. However, downstream of the bifurcation (B-B0 , centre, and C-C0 , right, cross-sections), the velocity profiles along H-H diameter are skewed towards the inner PRB wall due to the centrifugal forces exerted on the fluid. The deformation shifts further towards the inner wall of the PRB, in addition to the appearance of a second, much lower, peak

Table 2 Fluid and particle properties used in the PRB simulations. Fluid and particle properties Fluid temperature, T f Fluid density, qf Fluid dynamic viscosity,

293 K 1.21 kg/m3

lf

Particle density, qp

1:81  105 kg/m s 900 kg/m3

Particle diameter, dp

1–10 lm

Rest Reynolds number, Re Fluid mean inlet velocity, Fluid inlet flow, Q o Light exercise Reynolds number, Re Fluid mean inlet velocity, Fluid inlet flow, Q o Heavy exercise Reynolds number, Re Fluid mean inlet, velocity Fluid inlet flow, Q o

to

464 1.24 m/s 1.83 l/min

to

1132 3.02 m/s 4.46 l/min

to

1788 4.78 m/s 7.06 l/min

5598

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

Fig. 2. Definition of cross-sections and diameters of interest.

Fig. 3. Axial velocity profiles for different Reynolds numbers at cross-sections A-A (left), B-B0 (centre) and C-C0 (right) for symmetric flow along diameters: (a) H-H, and (b) V-V.

of the axial velocity towards the outer wall of the PRB with increasing Reynolds number. The axial velocity profiles along the diameters V-V for the A-A cross-section are still almost parabolic, whereas for the B-B0 and C-C0 cross-sections are deformed but remain symmetric with respect to the symmetry plane z ¼ 0 of the PRB. These results are consistent with other numerical solutions [40,41]. Likewise, in Fig. 4 the effect of flow asymmetry on axial velocity profiles is shown for both daughter tubes of the PRB along diameters H-H (top) and V-V (bottom) for Re ¼ 1132. In contrast to the symmetric flow cases, the axial velocity profile at the A-A cross-section along H-H diameter is no longer symmetric. The skewness of the curves is a result of the blockage effect in daughter tube 2 [41], and gets more pronounced in the totally obstructed flow case (Q 2 ¼ 0). Downstream the bifurcation, the axial velocity profiles along H-H are still shifted towards the inner PRB wall, but the flow in the second daughter tube is markedly less than the one in the first daughter tube. Moreover, it can be seen in the graphs along H-H diameter for

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

5599

Fig. 4. Axial velocity profiles for Re ¼ 1132, symmetric and asymmetric flow along diameters: (a) H-H, and (b) V-V. Corresponding cross-sections at both daughter tubes are shown.

the B1-B10 and C1-C10 that the fluid accelerates in daughter 1, comparing to the symmetric case, in order to accommodate the excess flow. The secondary-flow streamlines and contours of constant axial fluid velocity are shown in Fig. 5(a)–(c) at cross-sections B-B0 and C-C0 for the symmetric flow cases (left). Moreover, the contours of constant axial fluid velocity at the geometric symmetry plane z ¼ 0 and the velocity vectors at the H-H diameter of the cross-sections are depicted (right). The streamlines of the low Reynolds number secondary flow show the formation of a pair of symmetric, counter-rotating vortices. The centers of the vortices are slightly displaced towards the outer PRB wall at the C-C0 cross-section. Moreover, the peak of the axial fluid velocity is located closer to the inner wall. The secondary-flow streamlines for higher Reynolds numbers, Re ¼ 1132 and Re ¼ 1788, at C-C0 cross-section also show two main, symmetric counter-rotating vortices, but their centers are displaced towards the outer PRB wall and they are skewed with respect to the symmetry plane. In addition, increased centrifugal forces lead to increased fluid flow towards the inner wall. These results are quantitatively in agreement with other numerical simulations [41–44]. Similarly, the secondary flow streamlines and contours of constant axial fluid velocity downstream the bifurcation in both daughter tubes of the PRB are shown in Fig. 5 for the asymmetric flows Q 1 =Q 2 ¼ 2 (d) and Q 2 ¼ 0 (e). The flow asymmetry between the daughter tubes is obvious for both the fluid velocity vectors and the velocity magnitude contour. The maximum velocity at A-A cross-section is shifted towards daughter tube 1, the effect being more obvious for Q 2 ¼ 0. Comparing the corresponding cross-sections between the daughter tubes, one can notice the differences in the secondary flow, as well.

5600

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

Fig. 5. Contour of fluid velocity magnitude and vectors at the plane of symmetry of the PRB (z ¼ 0) and secondary flow streamlines and constant velocity magnitude contours at different cross-sections downstream the bifurcation: (a) Q 1 =Q 2 ¼ 1, Re ¼ 464, (b) Q 1 =Q 2 ¼ 1, Re ¼ 1132, (c) Q 1 =Q 2 ¼ 1, Re ¼ 1788, (d) Q 1 =Q 2 ¼ 2, Re ¼ 1132, and (e) Q 2 ¼ 0, Re ¼ 1132.

In cross-sections B1-B10 and C1-C10 of daughter 1 for both cases, there is the development of two counter-rotating vortices with their centers approaching the outer PRB wall, such as those for the symmetric case of Fig. 5(b), but the maximum axial fluid velocity is less shifted towards the inner PRB wall. These features also appear, less intense, in daughter tube 2 of the partially obstructed flow case, Q 1 =Q 2 ¼ 2, as shown in the right part of Fig. 5(d). However, this is not the case for daughter tube 2 of the totally obstructed flow, Q 2 ¼ 0. There are two counter-rotating vortices at the B2-B20 cross-section, but their centers are displaced toward the inner PRB wall, whereas there are no vortices at the C2-C20 cross-section. It should be noted that the white background color in these cross-sections (right part of Fig. 5(e)), indicates the almost zero fluid velocity in daughter tube 2 in this case. Closer observation of the fluid velocity vectors in Fig. 5, indicate that there are regions in the daughter tubes that flow recirculation occurs, e.g., at C–C0 location of Fig. 5(c) and and at the B2-B20 location Fig. 5(e). In Fig. 6, the streamlines of the main fluid flow at the plane of geometric symmetry, z ¼ 0, are shown for all cases. For the symmetric flow cases, Q 1 =Q 2 ¼ 1, there is recirculation for the higher Reynolds numbers, Re ¼ 1132 and Re ¼ 1788, and by optical comparison one notices that the recirculation zone becomes wider with increasing Reynolds number. On the other hand, for the asymmetric flows, recirculation occurs in the obstructed daughter tube, which for Q 2 ¼ 0 is very strong. 3.2. Aerosol particles The fraction of the deposited particles in the PRB, g, is calculated by:

_ cj

g ¼ 1  _ outlet ; cjinlet

ð16Þ

R where c_ ¼ St c~ v p  d~St is the dimensionless particle flowrate through the cross-sectioned area S~t . The subscripts ‘‘inlet’’ and ‘‘outlet’’ refer to the PRB inlet and outlet, respectively. The particle diameter ranges between 1 lm and 10 lm, for which inertial and gravitational effects should be taken into account in the calculations. The orientation of the PRB is such that the acceleration of gravity is equal to ~ g ¼ ðg x ; g y ; g z Þ ¼ ð0; 1; 0Þ, in order to coincide with the orientation of the experimental set-up of Kim and Iglesias [4], the findings of which are compared against the results of the present study. Total deposition fraction, g, for symmetric flow conditions (Q 1 =Q 2 ¼ 1) as a function of the particle Stokes number, St, is presented in Fig. 7 for different Reynolds numbers. The results of the Eulerian model are in good agreement with the exper-

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

5601

Fig. 6. Main fluid flow streamlines at the plane of symmetry of the PRB (z ¼ 0): (a) Q 1 =Q 2 ¼ 1, Re ¼ 464, (b) Q 1 =Q 2 ¼ 1, Re ¼ 1132, (c) Q 1 =Q 2 ¼ 1, Re ¼ 1788, (d) Q 1 =Q 2 ¼ 2, Re ¼ 1132, and (e) Q 2 ¼ 0, Re ¼ 1132.

Fig. 7. Total deposition fraction for symmetric flow. Comparison with experimental results [8].

imentally derived results [8] up to St ¼ 0:08 for the range of Reynolds number under study, whereas the present model underestimates deposition fraction for higher Stokes numbers. The effect of asymmetric fluid flow on the total deposition fraction is shown in Fig. 8. This results correspond to Reynolds number equal to 1132. From this figure it is shown that partially obstructed flow in daughter tube 2, which results in a Q 1 =Q 2 ¼ 2 flow ratio, does not affect significantly total deposition fraction in the PRB comparing with the symmetric flow

5602

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

Fig. 8. Total deposition fraction for symmetric and asymmetric flow (Re ¼ 1132).

case (Q 1 =Q 2 ¼ 1). This finding is consistent with the experimental results of Kim and Fisher [9], who studied flow asymmetry in a PRB cast of lung generations three to five (G3-G5). From Fig. 8, it is also obvious that if Q 2 ¼ 0 in the single PRB under study here, total deposition fraction is considerably lower than in Q 1 =Q 2 ¼ 1 or Q 1 =Q 2 ¼ 2 cases. The effect of flow Reynolds number under symmetric flow conditions on the concentration profiles and particle deposition sites are shown in Figs. 9(a)–(c) and 10(a)–(c) for low (St ¼ 0:01) and high (St ¼ 0:1) particle Stokes number, respectively. For St ¼ 0:01 (Fig. 9), deposition occurs mainly at the bifurcation and at the middle of the upper and lower PRB walls. Deposition is also high at the bifurcation in the case of high Stokes number (Fig. 10), but there is also significant

Fig. 9. Particle concentration profiles at cross-sections downstream the bifurcation and particle deposition sites on the PRB wall for low Stokes number, St ¼ 0:01: (a) Q 1 =Q 2 ¼ 1, Re ¼ 464, (b) Q 1 =Q 2 ¼ 1, Re ¼ 1132, (c) Q 1 =Q 2 ¼ 1, Re ¼ 1788, (d) Q 1 =Q 2 ¼ 2, Re ¼ 1132, and (e) Q 2 ¼ 0, Re ¼ 1132.

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

5603

Fig. 10. Particle concentration profiles at cross-sections downstream the bifurcation and particle deposition sites on the PRB wall for high Stokes number, St ¼ 0:1: (a) Q 1 =Q 2 ¼ 1, Re ¼ 464, (b) Q 1 =Q 2 ¼ 1, Re ¼ 1132, (c) Q 1 =Q 2 ¼ 1, Re ¼ 1788, (d) Q 1 =Q 2 ¼ 2, Re ¼ 1132, and (e) Q 2 ¼ 0, Re ¼ 1132.

deposition on the upper and lower PRB walls, which moves toward the outer wall and further downstream in the daughter tubes as Reynolds number increases. If the particle concentration profiles at the cross-sections B-B0 and C-C0 of these figures are observed in conjunction with the secondary flow streamlines of Fig. 6, it is shown that low Stokes number particles (St ¼ 0:01) accumulate at the periphery of the vortices, whereas higher Stokes number particles (St ¼ 0:1) are pushed further towards the PRB walls by the secondary flow. In fact, higher Reynolds number flows lead to accumulation of particles at the vicinity of the outer PRB wall, as the centers of the vortices in this case are also closer to the outer wall, and this is more pronounced further downstream the daughter tubes, i.e., at C-C0 cross-section rather than at B-B0 . Thus, the chance of a highly inertial particle to deposit at the outer wall of the PRB and towards the exit of the PRB increases monotonically with Reynolds number. Figs. 9(d) and (e) and 10(d) and (e), present the effect of flow asymmetry on particle deposition sites along the PRB and concentration profiles at cross-sections downstream the bifurcation for St ¼ 0:01 and St ¼ 0:1, respectively. The patterns for the partially obstructed flow (Q 1 =Q 2 ¼ 2), resemble those of the symmetric case with the same Reynolds number (Re ¼ 1132). However, deposition in daughter 2 occurs earlier along the tube and mostly on the upper and lower walls, in contrast with deposition in daughter 1, where peak deposition occurs towards the exit and the outer PRB wall. Moreover, the less intense secondary flow in daughter tube 2 (Fig. 6(d)) allows for the particles to accumulate at and outline better the periphery of the vortices, as it is shown by comparing cross-sections B1-B10 and C1-C10 with B2-B20 and C2-C20 , respectively, in Figs. 9(d) and 10(d). These features are particularly obvious for the higher Stokes number. In the case of Q 2 ¼ 0, it is shown that there is particle deposition in the totally obstructed daughter tube 2, which in the case of St ¼ 0:01 (Fig. 9(e)) is comparable to the deposition at the bifurcation. Particles that enter daughter 2 are pushed to the periphery of the recirculation zones, shown in Fig. 6(e), and probably some of them are trapped between them increasing the chances to deposit in the region. This is true for high Stokes number particles (Fig. 10(e)) as well, where deposition occurs also at the outer wall and towards the exit of daughter tube 1. Low Stokes number particles deposit in a wider region of the PRB daughter tubes including the inner wall (Fig. 9), whereas there are regions on the outer wall that are almost free of particles. On the other hand, deposition for high Stokes particles (Fig. 10) is more localized towards the outer wall and less particles deposit on the inner PRB wall. In all studied cases, there is a region of the outer wall, at the beginning of the daughter tubes opposite the bifurcation, where almost no particle deposition occurs. Finally, it should be noticed that even if the total deposition fraction does not differs greatly between the various flow cases, with the Q 2 ¼ 0 case exception, the particle deposition sites and concentration profiles are affected by the flow conditions.

5604

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

4. Conclusions A computational study of particle transport and deposition in a physiologically realistic bifurcation corresponding to generations three and four (G3-G4) of the human lung was conducted using a fully Eulerian model. In particular, the effect of particle size and fluid flow Reynolds number and symmetry, was investigated. It was found that total particle deposition fraction is loosely dependent on fluid flow Reynolds number and our results are in agreement with the experimental findings [8]. Total deposition fraction, also, does not change significantly between symmetric Q 1 =Q 2 ¼ 1 and asymmetric Q 1 =Q 2 ¼ 2 flow conditions, but is considerably lower for the totally obstructed case, Q 2 ¼ 0. Even though deposition fraction did not vary much for the various flow conditions, deposition sites and particle concentration profiles depend on them. Apart from deposition at the bifurcation, which was present at all cases, deposition sites downstream the bifurcation move towards the outer PRB wall and daughter tube exit as Reynolds number increases. For the Q 2 ¼ 0, it was also shown that particles are trapped and deposit even in the obstructed daughter tube. In all studied cases, it was found that there is a particle-free region at the outer wall, at the beginning of the daughter tubes opposite the bifurcation. In the present work, a single bifurcation is considered, at the inlet of which a parabolic air-flow profile is assumed. This assumption permitted the comparison with the measurements of Kim et al. [8]. However, this is not the case in reality, as the upstream bifurcations would result in different inlet flow profiles. The more realistic description of the air flow field in the PRB taking into account the effect of previous generations it was beyond the scope of this study and should be considered seperately. The use of the present Eulerian model, which is a CFD-based fully mechanistic model, to obtain not only the deposition fraction, but also the characteristics of particle deposition and concentration patterns, can offer useful physical insight under different aerosol flow conditions, in a straightforward and simple manner. At the same time, however, it is rather complicated and requires considerable computing power. As a result, it is applicable to selected structural elements (generations) of the respiratory tract. On the other hand, there are simple, validated, operational models, which use a combination of theoretical and empirical expressions to predict particle deposition in the whole respiratory tract. One such model is MPPD (Multiple Path Particle Dosimetry) model [45], a computational model for calculating the deposition and clearance of monodispersed and polydispersed aerosols in the respiratory tract of rats and humans (adults and children). The model deals with particles sizes ranging from 0.01 lm to 20 lm, and is based both on single- and multiple-path methodologies for the determination of air flow and aerosol deposition in the lung. The particles deposition mechanisms taken into account in the model are diffusion, sedimentation and impaction, and within each airway deposition is calculated using theoretical deposition efficiencies for each one of these mechanisms. Moreover, empirical efficiency functions are used to determine aerosol filtration by the head, whereas the three-compartment model of ICRP (1994) and a two-compartment model with semi-empirical clearance rate constants are used to calculate clearance for humans and rats, respectively. As all semi-empirical models, however, it is limited to the specific morphology, physiology and lung conditions for which its parameters were obtained, and also neglects aerosol dynamics [46]. Nevertheless, one can combine the positive features of each model; a well validated semi-empirical model of the whole respiratory tract, such as MPPD, can be used as screening method in order to recognize the generations of the lung that deposition peaks, and then a detailed three-dimensional CFD model can be employed to study both the air flow field and particle transport and deposition in the generations of interest in full detail. Acknowledgement This study was partially supported by project NanoTEST under Contract No. HEALTH-2008-201335 of the European Commission. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

A. Annapragada, N. Mishchiy, In silico modeling of aerosol deposition in lungs, Drug Discov. Today Dis. Models 4 (2007) 155–161. R.B. Schlesinger, D.E. Bohning, T.L. Chan, Particle deposition in a hollow cast of the human trancheobronchial tree, J. Aerosol Sci. 8 (1977) 429–445. Y.-S. Cheng, Y. Zhou, B.T. Chen, Particle deposition in a cast of human oral airways, Aerosol Sci. Technol. 31 (1999) 286–300. C. Kim, A. Iglesias, Deposition of inhaled particles in bifurcating airway models with varying airway geometry, J. Aerosol Med. 2 (1989) 1–14. D.Y.H. Pui, F. Romay-Novas, B.Y.H. Liu, Experimental study of particle deposition in bends of circular cross section, Aerosol Sci. Technol. 7 (1987) 300– 315. F. Cai, C. Yu, Inertial and interceptional deposition of spherical particles and fibers in a bifurcating airway, J. Aerosol Sci. 19 (1988) 679–688. Y.S. Cheng, C.S. Wang, Motion of particles in bends of circular pipes, Atmos. Environ. 15 (1981) 301–306. C.S. Kim, D.M. Fisher, D.J. Lutz, T.R. Gerrity, Particle deposition in bifurcating airway models with varying airway geometry, J. Aerosol Sci. 25 (1994) 567–581. C.S. Kim, D.M. Fisher, Deposition characteristics of aerosol particles in sequentially bifurcating airway models, Aerosol Sci. Technol. 31 (1999) 198–220. M.J. Oldham, R.F. Phalen, T. Heistracher, Computational fluid dynamic predictions and experimental results for particle deposition in an airway model, Aerosol Sci. Technol. 32 (2000) 61–71. A.F. Miguel, A.H. Reis, M. Aydin, Aerosol particle deposition and distribution in bifurcating ventilation ducts, J. Hazardous Mater. 116 (2004) 249–255. J. Lee, J. Goo, Numerical simulation of airflow and inertial deposition of particles in a bifurcation channel of square cross-section, J. Aerosol Med. 5 (1992) 131–154. B. Asgharian, S. Anjilvel, Inertial and gravitational deposition of particles in a square cross section bifurcating airway, Aerosol Sci. Technol. 20 (1994) 177–193.

M. Pilou et al. / Applied Mathematical Modelling 37 (2013) 5591–5605

5605

[14] W. Hofmann, R. Golser, I. Balásházy, Inspiratory deposition efficiency of ultrafine particles in a human airway bifurcation model, Aerosol Sci. Technol. 37 (2003) 988–994. [15] J. Comer, C. Kleinstreuer, C. Kim, Flow structures and particle deposition patterns in double-bifurcation airway models. Part 2. Aerosol transport and deposition, J. Fluid Mech. 435 (2001) 55–80. [16] P.W. Longest, M.J. Oldham, Mutual enhancements of CFD modeling and experimental data: A case study of 1 lm particle deposition in a branching airway model, Inhal. Toxicol. 18 (2006) 761–771. [17] P.W. Longest, S. Vinchurkar, Inertial deposition of aerosols in bifurcating models during steady expiratory flow, J. Aerosol Sci. 40 (2009) 370–378. [18] P.W. Longest, J. Xi, Computational investigation of particle inertia effects on submicron aerosol deposition in the respiratory tract, J. Aerosol Sci. 38 (2007) 111–130. [19] Z. Zhang, C. Kleinstreuer, J.F. Donohue, C.S. Kim, Comparison of micro- and nano-size particle depositions in a human upper airway model, J. Aerosol Sci. 36 (2005) 211–233. [20] Z. Zhang, C. Kleinstreuer, C.S. Kim, Gas-solid two-phase flow in a triple bifurcation lung airway model, Int. J. Multiphase Flow 28 (2002) 1021–1046. [21] A. Farkas, I. Balásházy, Quantification of particle deposition in asymmetrical tracheobronchial model geometry, Comput. Biol. Med. 38 (2008) 508–518. [22] Z. Zhang, C. Kleinstreuer, C.S. Kim, Comparison of analytical and cfd models with regard to micron particle deposition in a human 16-generation tracheobronchial airway model, J. Aerosol Sci. 40 (2009) 16–28. [23] G. Yu, Z. Zhang, R. Lessmann, Computer simulation of the flow field and particle deposition by diffusion in a 3-D human airway bifurcation, Aerosol Sci. Technol. 25 (1996) 338–352. [24] G. Yu, Z. Zhang, R. Lessmann, Fluid flow and particle diffusion in the human upper respiratory system, Aerosol Sci. Technol. 28 (1998) 146–158. [25] Z. Zhang, C. Kleinstreuer, Species heat and mass transfer in a human upper airway model, Int. J. Heat Mass Tranfer 46 (2003) 4755–4768. [26] P.W. Longest, M.J. Oldham, Numerical and experimental deposition of fine respiratory aerosols: Development of a two-phase drift flux model with near-wall velocity corrections, J. Aerosol Sci. 39 (2008) 48–70. [27] L. Gradon, D. Orlicki, Deposition of inhaled aerosol particles in a generation of the tracheobronchial tree, J. Aerosol Sci. 21 (1990) 3–19. [28] E.R. Weibel, Morphometry of the Human Lung, Springer-Verlag, 1963. [29] M. Pilou, S. Tsangaris, P. Neofytou, C. Housiadas, Y. Drossinos, Inertial particle deposition in a 90-degree laminal-flow bend: An Eulerian fluid-particle approach, Aerosol Sci. Technol. 45 (2011) 1376–1387. [30] P. Neofytou, S. Tsangaris, Flow effects of blood constitutive equations in 3D models of vascular anomalies, Int. J. Numer. Meth. Fluids 51 (2006) 489– 510. [31] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, Taylor & Francis, 1980. [32] J. Fernandez de la Mora, D.E. Rosner, Effects of inertia on the diffusional deposition of small particles to spheres and cylinders at low Reynolds numbers, J. Fluid Mech. 125 (1982) 379–395. [33] Y. Drossinos, C. Housiadas, Aerosol flows, in: C.T. Crowe (Ed.), Multiphase Flow Handbook, CRC Press, Boca Raton, 2006. [34] J.H. Ferziger, M, edition3rd edition., Peric´, Computational Methods for Fluid Dynamics, Springer, 2002. [35] T. Heistracher, W. Hofmann, Physiologically realistic models of bronchial airway bifurcations, J. Aerosol Sci. 26 (1995) 497–509. [36] E. Makris, V. Gkanis, S. Tsangaris, C. Housiadas, A methodology to generate structured computational grids from DICOM data: application to a patientspecific abdominal aortic aneurysm (AAA) model, Comput. Methods Biomech. Biomed. Engin., vol. 15, 2011, p. 173. [37] E. Makris, S. Tsangaris, C. Housiadas, A methodology to generate a patient specific high quality structured computational domain from medical imaging data, in: Proceedings of the ECCOMAS Thematic International Conference on Simulation and Modeling of Biological Flows (SIMBIO 2011). [38] P.W. Longest, S. Vinchurkar, Validating CFD predictions of respiratory aerosol deposition: Effects of upstream transition and turbulence, J. Biomech. 40 (2007) 305–316. [39] E. Makris, P. Neofytou, S. Tsangaris, C. Housiadas, A novel method for the generation of multi-block computational structured grids from medical imaging of arterial bifurcations, Inpress, Med. Engin. Phys, 2011. [40] I. Balásházy, W. Hofmann, Particle deposition in airway bifurcations - I. Inspiratory flow, J. Aerosol Sci. 24 (1993) 745–772. [41] J. Comer, C. Kleinstreuer, Z. Zhang, Flow structures and particle deposition patterns in double-bifurcation airway models. Part 1. Air flow fields, J. Fluid Mech. 435 (2001) 25–54. [42] P.W. Longest, S. Vinchurkar, Effects of mesh style and grid convergence on particle deposition in bifurcating airway models with comparisons to experimental data, Med Eng. Phys. 29 (2007) 350–366. [43] T.B. Martonen, X. Guan, R.M. Schreck, Fluid dynamics in airway bifurcations: II. Secondary currents, Inhal. Toxicol. 13 (2001) 281–289. [44] T.B. Martonen, X. Guan, R.M. Schreck, Fluid dynamics in airway bifurcations: III. Localized flow conditions, Inhal. Toxicol. 13 (2001) 291–305. [45] O. Price, B. Asgharian, F. Miller, F. Cassee, R. de Winter-Sorkina, Multiple Path Particle Dosimetry model (MPPD v1.0), CIIT & RIVM, 2002. [46] C. Housiadas, M. Lazarides, Inhalation dosimetry modelling, in: M. Lazarides, I. Colbeck (Eds.), Human Exposure to Pollutants via Dermal Absorption and Inhalation, Springer, Netherlands, 2010.