Large eddy and detached eddy simulations of fluid flow and particle deposition in a human mouth–throat

Aerosol Science 39 (2008) 862 – 875 www.elsevier.com/locate/jaerosci

Large eddy and detached eddy simulations of fluid flow and particle deposition in a human mouth–throat S.T. Jayarajua,∗ , M. Brounsa , C. Lacora , B. Belkassemb , S. Verbanckc a Department of Mechanical Engineering, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussel, Belgium b Department of Mechanics of Materials and Constructions, Vrije Universiteit Brussel, Pleinlaan 2, 1050 Brussel, Belgium c Respiratory Division, Academic Hospital, Vrije Universiteit Brussel, Laarbeeklaan 101, 1090 Brussel, Belgium

Received 22 February 2008; received in revised form 29 May 2008; accepted 2 June 2008

Abstract Fluid flow was simulated in a human mouth–throat model under normal breathing conditions (30 l/min) alternatively employing RANS k . (without near-wall corrections), DES and LES methods. To test the validity of the fluid phase simulations, PIV measurements were carried out in a central sagittal plane of an identical model cast. Velocity and kinetic-energy profiles showed good quantitative agreement of experiments with LES/DES, and less so with RANS k .. Mouth–throat deposition was simulated for particle diameters 2, 4, 6, 8 and 10 m. By comparison with existing experimental data, LES/DES showed considerable improvement over the RANS k . model in predicting deposition for particle sizes below 5 m. For the bigger particles, RANS k . and LES/DES methods produced similarly good predictions. It is concluded that for the simulation of medication aerosols inhaled at a steady flow rate of 30 l/min, LES and DES provide more accurate results than the RANS k . model tested. 䉷 2008 Elsevier Ltd. All rights reserved. Keywords: Mouth–throat model; Particle image velocimetry (PIV); Large eddy simulation (LES); Detached eddy simulation (DES); Particle deposition

1. Introduction Inhaled medication is a preferred method of drug administration to the lung for the first-line therapy of asthma and chronic obstructive pulmonary diseases. The inhaled aerosol particles need to negotiate the mouth–throat structure in order to reach the smaller airways and the alveolar lung zone that could benefit from aerosol therapy. The complexity of the extrathoracic portion of the oral airway, which includes bends and sudden cross-sectional changes, potentially induces considerable local medication deposition before actually reaching the lungs. A quantitative study of aerosol transport and deposition in a realistic model of the mouth–throat poses many challenges, both from an experimental and simulation standpoint. Most previous experimental studies in mouth–throat geometries have focussed on the characterization of fluid flow fields, with the ultimate aim of predicting aerosol deposition patterns. For instance, Corcoran and Chigier (2000) used phase doppler interferometry to study the axial flow field in a larynx/trachea cast, assessing the effects of the laryngeal ∗ Corresponding author. Tel.: +32 2 6292336; fax: +32 2 6292880.

E-mail addresses: [email protected], [email protected] (S.T. Jayaraju). 0021-8502/$ - see front matter 䉷 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.jaerosci.2008.06.002

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

863

jet on inhalation flow patterns. Gemci, Corcoran, and Chigier (2002) used laser doppler velocimetry to characterize axial velocity fields and turbulent intensity levels in a simple throat model (essentially a constricted tube). Heenan, Matida, Pollard, and Finlay (2003) used endoscopic particle image velocimetry (PIV) to visualize the velocity fields in distinct parts of the central sagittal plane of an idealized model of the human oropharynx. In specific locations of the same geometry, Johnstone, Heenan, Uddin, Pollard, and Finlay (2004) measured mean and RMS axial velocity using x-hot-wire anemometry. Heenan, Finlay, Grgic, Pollard, and Burnell (2004) combined PIV and scintigraphic aerosol deposition measurements in an attempt to relate time-averaged flow fields in the central sagittal plane to local aerosol deposition in two realistic extrathoracic airway geometries. Finally, DeHaan and Finlay (2001) employed ultraviolet spectrophotometric assay to measure depositions in a mouth–throat geometry for different inhalation devices, and DeHaan and Finlay (2004) collected aerosols from dry powder inhalers on disposable filters to determine oral cavity deposition. Most previous numerical studies have simulated the fluid/particle behavior in mouth–throat geometries by means of RANS simulations, mainly using two equation turbulence models like k. and k. for the fluid phase, sometimes coupled with eddy interaction model (EIM) for the particle phase. Stapleton, Guentsch, Hoskinson, and Finlay (2000) observed that even in a simplified mouth–throat model, k. turbulence model was not suitable for the accurate prediction of particle deposition. Matida, Finlay, Lange, and Grgic (2004) obtained better results with standard k. model compared to standard k., yet, for the particles in low Stokes number range (pertinent to medical aerosols), the simulated total mouth–throat deposition was ∼ 50% in contrast to the much lower depositions encountered experimentally (i.e., less than 5%). When anisotropy effects close to the walls were taken into account in the EIM, simulated deposition came down to ∼ 15%. On the other hand, Xi and Longest (2007) employed a low Reynolds number k. model to assess the effects of geometry simplifications on aerosol deposition, reporting deposition values of less than 20% (even without considering anisotropy near the walls). RANS models, which have been basically developed for fully turbulent flows, may be inappropriate for the prediction of particle deposition in the mouth–throat, where flow is in fact transitional. Realizing this problem, several authors have started to explore the possibility of using large eddy simulation (LES) methods for the study of particle deposition in the mouth–throat, by applying LES in relatively simple structures such as turbulent duct flow (Tian & Ahmadi, 2007), a 90◦ bend (Breuer, Baytekin, & Matida, 2006) or a constricted tube (Luo, Hinton, Liew, & Tan, 2004). LES of deposition in a simplified mouth cavity model (Breuer, Matida, & Delgado, 2007; Ilie, Matida, & Finlay, 2008; Matida, Finlay, Breuer, & Lange, 2006) and in a simplified mouth–throat model (Jin, Fan, Zeng, & Cen, 2007) have indicated the potential of LES to more accurately simulate aerosol transport in the extrathoracic airways. The aim of the present work was to test the validity of the most commonly used modeling methods, namely the RANS, LES and detached eddy simulation (DES) for the description of fluid/particle behavior in the mouth–throat model. The existing experimental data either resulted from intrusive measurement methods or only provided partial data sets on existing model geometries making the comparison with present simulations difficult. Therefore, we perform PIV measurements in a three-dimensional mouth–throat cast and use exactly the same mouth–throat geometry for CFD simulations, enabling direct comparison between simulations and experiments. On the part of CFD simulations, RANS employed SST k. turbulence model, DES was based on the Spalart–Allmaras model for the near-wall region, and LES was based on two sub-grid-scale models, namely the Smagorinsky–Lilly and the WALE model. Alternatively, we considered the frozen LES method proposed by Matida et al. (2006) where the particles are released in a frozen (static) instantaneous velocity field and tracked in steady mode, without any EIM. 2. Experimental methods The mouth–throat geometry used for PIV measurements and CFD simulations is shown in Fig. 1. This geometry was previously used for the CFD study of tracheal stenosis (Brouns et al., 2007). From this geometry, a male model was generated by means of stereolithography, after which a one block transparent female model of the mouth– throat was obtained by following the procedures which are discussed in detail by Hopkins, Kelly, Wexler, and Prasad (2000). No geometric re-scaling was necessary for the present PIV measurements. A water–glycerin mixture was used with a viscosity of 5.88 × 10−6 m2 /s as determined with the AVS300 viscosimeter from Schott Gerate at 25.2 ◦ C. Dynamic similarity was achieved for the liquid-based experiments by matching the Reynolds number. The PIV measurements were performed in a central sagittal plane at a volumetric flow rate of 11.2 l/min, corresponding to a Re number of 2527 (based on inlet diameter).

864

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

Fig. 1. Simplified geometry reconstructed from CT-scan data showing different cross-sections of the geometry.

The water–glycerin mixture was pumped in a reservoir placed ∼ 1.5 m above the model (Fig. 2) in order to create a developed velocity profile. The reservoir had an inlet which was connected to the pump, an outlet connecting the model and an overflow exit which guaranteed a constant level in the reservoir. The outlet of the reservoir was separated from the inlet and overflow exit by a fine maze to stabilize the level and to remove any fluctuations caused by the pump. The flow rate was regulated by a valve which was placed behind the flow meter (KOBOLD Instrumentation NV/SA with accuracy of 4% f.s.). A new wave MinilaseII Nd-Yag laser (532 nm wavelength, 100 mJ/pulse) was synchronized with a pulse separation, depending on the flow rate. The pulse separation was chosen in such a way that the reflection of the tracer particles (10 m hollow glass spheres) shift 5 pixels between an image pair. The laser beams were combined and formed into a sheet with cylindrical optics. This pulsed sheet was passed through the model, parallel to the flow, and the light scattered from the particles was recorded with a PCO sensicam QE 5 Hz camera.

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

865

Fig. 2. Schematic representation of the experimental setup.

Approximately 4000 image pairs were recorded. The images were analyzed using PIVview 2C software (PIVTEC GmbH, Germany). The vector fields were generated using cross-correlation fast Fourier transform (FFT) with a multigrid procedure combined with a sub-pixel-based image shifting or image deformation with a third-order interpolation scheme (Scarano & Riethmuller, 2000). The final interrogation region was 32 × 32 pixels with an overlap of 50%. Spurious vectors were detected by using the normalized median test, which eliminates a dependence of the detection criterion on the interrogation domain size (Westerweel & Scarano, 2005). Only 0.1% of the vectors had to be interpolated. A least squares 3-point Gauss fit algorithm was used to recover the sub-pixel displacement of the correlation. 3. Numerical methods The fluid and particle phase were solved employing the incompressible solver of FLUENT 6.3. 3.1. Fluid phase RANS: The time-averaged Navier–Stokes equations are modeled employing low Reynolds number variant of SST k. model (Menter, 1994), which requires resolving the near-wall region with a fine mesh. This model has been selected

866

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

based on its ability to accurately predict the particle depositions in the models of mouth–throat geometries (Kleinstreuer & Zhang, 2003; Matida et al., 2004; Xi & Longest, 2007). Second-order upwind scheme for momentum equation and third-order MUSCL scheme for k. equation were employed for spatial discretization. SIMPLE algorithm was used for pressure–velocity coupling. DES: DES is most often referred to as the hybrid RANS/LES model where unsteady RANS is employed to model the near-wall region and LES in the core turbulent region. The present DES model is based on the standard one equation Spalart–Allmaras model. The length scale in the DES model is defined as d˜ = min(d, Cdes ), where d is the distance to closest wall and  is the grid spacing based on the largest grid space in the x, y or z directions forming the computational cell. The empirical constant Cdes has a value of 0.65. Second-order implicit formulation for temporal discretization and central differencing for spatial discretization of momentum as well as the turbulent viscosity equation were employed. DES is computationally more expensive than RANS but less expensive than LES since the near wall is modeled using RANS approach. However, DES involves solving of an additional turbulent viscosity equation. LES: In LES, the big three-dimensional eddies which are dictated by the geometry and boundary conditions of the flow involved are directly resolved, whereas the small eddies which tend to be more isotropic and less dependent on the geometry are modeled. The sub-grid-scale stresses (ij ) resulting from the filtering operation is given by ij − 13 kk ij = −2t S ij

(1)

t is the sub-grid-scale turbulent viscosity. S ij is the rate of strain tensor. Two constant sub-grid-scale models, namely the Smagorinsky–Lilly model and the wall adapting local eddy viscosity (WALE) model were tested. In the Smagorinsky–Lilly model, the eddy-viscosity is modeled by  t = L2s 2S ij S ij (2) Ls is the mixing length for sub-grid scale given by min( d, Cs V 1/3 ). is the von Karman constant, d is the distance to the closest wall and Cs is the Smagorinsky constant which is taken to be 0.1. V is the computational cell volume. In the WALE model, the eddy viscosity is modeled by t = Ł2s

(Sijd Sijd )3/2 (S ij S ij )5/2 + (Sijd Sijd )5/4

(3)

Ls is given by min( d, C V 1/3 ). C has a default value of 0.325. The constant Cs in the Smagorinsky–Lilly model can also be dynamically computed based on the information provided by the resolved scales of motion. However, the influence of the sub-grid-scale model is expected to be small due to the low Reynolds number of flows we are dealing with in the present investigation (Breuer, 1998, 2000). Similar to DES, second-order implicit formulation is used for temporal discretization and central differencing for spatial discretization of momentum equation. 3.2. Particle phase Assuming large particle-to-air density ratio, negligible particle rotation, no inter-particle collision, and drag force as the dominant point force, the Lagrangian equations governing the particle motion are given by dxp = up dt

(4)

gx (p − ) dup = Fd (u − up ) + dt p

(5)

xp is position of the particle.  and p are the density of fluid and the particle, respectively. gx is the gravitational force which is oriented in the vertical direction. It should be noted that, for numerical simulations, the mouth–throat geometry is turned by 90 ◦ with respect to our experimental setup, in order to allow direct comparison with the experimental curve-fit of Grgic, Finlay, Burnell, and Heenan (2004).

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

867

The unbalanced pressure distribution on the surface of the particle due to the difference in fluid velocity (u) and particle velocity (up ) is termed as the drag force (Fd ) which is given by Fd =

18 Cd Re p dp2 24

(6)

 is the dynamic viscosity of the fluid. dp is the diameter of the particle. Cd is the coefficient of drag given by Cd = a1 +

a3 a2 + 2 Re Re

(7)

The ai coefficients are constants for smooth spherical particles as given by Morsi and Alexander (1972). Re is the particle Reynolds number defined as dp |up − u|/. RANS: The instantaneous velocity (u) comprises of a mean component (u) and a fluctuating component (u ). Assuming isotropic turbulence, the fluctuating component is expressed as (Gosman & Ioannides, 1981)  2  u = k×

(8) 3 where k is the turbulent kinetic energy of the flow and is a random number drawn from a Gaussian probability density function with zero mean and unit standard deviation. The chosen fluctuation is referred to a turbulent eddy whose time scale is given by Te = 2TL TL = CL

1 

(9) (10)

TL is the Lagrangian time scale and the constant CL for the SST k. turbulence model is 0.15. It is accepted that the fluctuation remains constant within a turbulent eddy during its life-time (Te ), while the respective mean velocity component (u) is varied according to particle position. The crossing trajectory effect is taken into account by limiting the crossing time to    Le tcross = − ln 1 − (11) |u − up |  is the relaxation time of the particle. DES and LES: In the unsteady mode, each fluid phase iteration is followed by a particle phase iteration and the particles are tracked in the real time. Hence, the effect of resolved large-scale instantaneous velocity on the particles are accounted for and there is no need for an EIM as in RANS. In case of frozen LES, the fluid flow is simulated for certain number of through flow cycles followed by the injection of particles in a stagnant (frozen) instantaneous velocity field and tracked as in RANS, but without the EIM. 4. Quality control The air flow in the mouth–throat geometry is investigated at a normal breathing rate of 30 l/min. A steady top hat velocity profile along with 5% turbulence intensity at the inlet and static pressure at the outlet were imposed. No-slip boundary condition was used at the walls. We investigated Lagrangian particles with density p = 912 kg/m3 and diameters 2, 4, 6, 8 and 10 m. Particular care was taken to obtain a uniform surface area distribution of particles at the model inlet. The particles were injected with their initial velocity set equal to that of inlet fluid velocity. As we are dealing with dilute particulate flow, one-way coupling as described by Elghobashi (1994) is assumed. Since the airway passage is normally wet, it is assumed that a particle is taken to be deposited on the wall as soon as it touches the wall. The effect on inlet condition on the fluid flow and particle deposition was investigated. Applying either a blunt inlet profile or a parabolic inlet profile showed negligible effect on the flow while comparing different cross-sectional velocities in the central sagittal plane of the model. A total deposition variation of less than 2% was observed for all particle diameters considered in the present study.

868

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

2 1

u/uin

u/uin

1.5

Experiment LES − 1.9*106 LES − 2.9*106 DES − 1.2*106

0.5

1 0.5

0

0 0

0.5 x/d

1

0

0.2

0.4

0.6

0.8

1

0.6

0.8

1

x/d

1.5

2 1.5 u/uin

u/uin

1 1

0.5 0.5 0

0 0

0.2

0.4

0.6 x/d

0.8

1

0

0.2

0.4 x/d

Fig. 3. Comparison of normalized 2-component (ux and uz ) velocity magnitude corresponding to the central sagittal plane. (a) Five millimeters above epiglottis (corresponds to section C in Fig. 1). (b–d) One, two and three tracheal diameters downstream of larynx, respectively (corresponds to section H, I, J in Fig. 1); x/d = 1 corresponds to the anterior airway wall.

RANS: A recent study by Vinchurkar and Longest (2008) has shown clear advantages of using hexahedral meshes in respiratory aerosol transport and deposition. The present mouth–throat model was meshed into 800,000 hexahedral cells with clustering in the vicinity of the wall and a stretching ratio of 1.2. The y + of the first layer of cells next to the wall was about 1. For each particle diameter, 5000 particles were injected at the model inlet. Meshes of either 800,000 or 2,170,000 cells (essentially obtained by doubling the number of grid points in each direction) had shown negligible variation of cross-sectional velocity magnitude in the larynx region and a deposition variation of less than 2% for all particle diameters. The difference in percentage deposition by injecting 5000 or 15,000 particles at the model inlet had been found to be less than 1%. Matida et al. (2004) had previously reported that their simulated deposition results were unaffected by increasing the number of injected particles from 1000 to 10,000. For the particle transport, an efficient automated tracking scheme, which is a combination of high order trapezoidal scheme and low order implicit scheme was employed. The typical simulation time of the flow field to obtain a convergence level of three orders of magnitude and to subsequently track 5000 particles took ∼ 14 h on an AMD Opteron 2.4 MHz dual-core processor. LES and DES: In case of LES, the mouth–throat model was meshed into 1.9 × 106 hexahedral cells with clustering towards the wall. The first cell layer next to the wall had a y + ∼ 0.2 with a stretching ratio of 1.05, which was sufficient to resolve the viscous sublayer. For the DES mesh, 1.2 × 106 hexahedral cells were employed with a y + ∼ 1 and a stretching ratio of 1.2. The converged steady state solution based on SST k. model was perturbed by adding random fluctuations and was used as an initial solution for faster convergence of LES/DES computations. To get rid of any possible initial condition effects, three through flow cycles were performed before starting the time-averaging. A through flow cycle is defined as the ratio between length of the airway model to the average cross-sectional velocity in the central sagittal plane. For 30 l/min, it is ∼ 0.1 s. A through flow cycle is much smaller compared to a typical inhalation period which is ∼ 2 s

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

869

(Li, Kleinstreuer, & Zhang, 2007). The time step to advance the flow was chosen such that the CFL number in the entire domain was less than 1, which ensures time advancement stability and accuracy. The typical time step was 1 × 10−5 s when using 1.9 × 106 mesh points. To test mesh independence of the solution, an additional simulation was performed for LES (WALE model) with the same time step but on 2.9 × 106 cells. Essentially, no differences were found in the average velocity magnitude except for a slight offset in the location of velocity profile on the anterior side at 5 mm above epiglottis (Fig. 3(a)). For all LES and DES simulations, obtaining a time-independent averaged solution took approximately four through flow cycles. Each through flow cycle for LES and DES took approximately 160 and 110 h, respectively, while running parallel on four AMD Opteron 2.4 MHz dual-core processors. Once the time-independent average solution was obtained, Lagrangian particles were introduced into the domain by injecting 5000 particles for each particle diameter uniformly over a time span of one through flow cycle. Previous LES studies (Jin et al., 2007; Matida et al., 2006) had indicated that ∼ 3500 particles were sufficient. Approximately five through flow cycles were required to get all the particles to either deposit on the wall or reach the outlet. Once the unsteady tracking was finished, the flow solution was also used for the frozen LES where the particles were injected in the frozen instantaneous velocity field and tracked as in RANS, but without the EIM. 5. Results and discussion 5.1. Fluid phase The flow patterns obtained with LES, DES and RANS are assessed by comparing the 2-component normalized velocity magnitude profiles with experimental ones at various model cross-sections (Fig. 4). Both LES sub-grid-scale

2

2 Experiment LES − Smagorinsky LES − Wale DES RANS

1.5 u/uin

u/uin

1.5

1

0.5

0.5

0

0 0

0.5 x/d

1

2

2

1.5

1.5 u/uin

u/uin

1

1

0.5

0

0.5 x/d

1

0

0.5 x/d

1

1

0.5

0

0 0

0.5 x/d

1

Fig. 4. Comparison of normalized 2-component (ux and uz ) velocity magnitude corresponding to the central sagittal plane. (a) Five millimeters above epiglottis. (b–d) One, two and three tracheal diameters downstream of larynx, respectively.

870

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

Experiment LES − Smag LES − Wale DES

0.2 0.1

0 0.5 x/d

1

0.3

0.3

0.2

0.2

2 k/uin

2 k/uin

0

0.1

0

0.5 x/d

1

0

0.5 x/d

1

0

0.5 x/d

1

0

0.5 x/d

1

0.1

0

0 0

0.5 x/d

1

0.3

0.3

0.2

0.2

2 k/uin

2 k/uin

0.2 0.1

0

0.1

0.1

0

0 0

0.5 x/d

1

0.3

0.3

0.2

0.2

2 k/uin

2 k/uin

LES − Smag (3 comp) RANS

0.3 2 k/uin

2 k/uin

0.3

0.1

0.1

0

0 0

0.5 x/d

1

2 Fig. 5. Left: comparison of normalized 2-component (u2 x and uz ) kinetic energy corresponding to the central sagittal plane between experiments, LES and DES. Right: comparison of normalized 3-component kinetic energy between LES and RANS. (a) Five millimeters above epiglottis. (b–d) One, two and three tracheal diameters downstream of larynx, respectively.

models (Smagorinsky and WALE) perform well in all four cross-sections, with a slightly better prediction of the velocity profile by the Smagorinsky model in the pharynx region on the anterior side (close to x/d = 1 in Fig. 4(a)). In all four cross-sections, DES and LES predictions of velocity profiles are very similar. By contrast, the widely used k. RANS model systematically overestimates velocity near the anterior airway wall, particularly in the tracheal region (Fig. 4(b–d)). In order to compare kinetic energy across the computations and experimental data, we first considered the 2-component kinetic-energy as measured by PIV. The experimental data were compared with the corresponding

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

871

Fig. 6. Left: time-averaged central sagittal plane velocity magnitude and corresponding streamlines. Right: time-averaged velocity magnitude (above 50% of maximum velocity in the airway model) and corresponding secondary velocity vector lines at six different cross-sections.

2-component kinetic-energy obtained with LES and DES (left panels of Fig. 5), showing a good agreement between simulations and experiments and negligible variation across the LES/DES models. Since the kinetic energy in RANS implicitly comprises of all three components, direct comparison with experimentally measured 2-component kinetic en2 ergy (u2 x and uz ) is not possible. Considering the similarity of 2-component kinetic-energy profiles between LES/DES and the experiments, the corresponding 3-component kinetic-energy profiles for LES Smagorinsky was compared with RANS (right panels in Fig. 5). It can be observed that the kinetic-energy profile obtained with the k. model is very different from that of LES, with both under- and overestimation of local kinetic energy depending on the cross-section under study. The time-averaged flow patterns and the secondary flow structures are presented in Fig. 6 for the LES Smagorinsky model. The six cross-sectional views show velocity magnitudes that are greater than 50% of the maximum velocity in the airway model (i.e., 2.6 m/s and above). The flow entering through the mouth piece impinges on the tongue and takes a bend upwards. As it continues to move forward, it accelerates in the middle part of the mouth due to reduction in cross-sectional area. As can be seen in sections A1–A2, velocity in most of this cross-section is above 2.6 m/s. Also, two distinct recirculation zones are seen. Towards the end of mouth region, the flow takes a downward turn and enters the pharynx in the form of a jet which undergoes an expansion due to increase in cross-sectional area. Consequently, the velocity is reduced and complex secondary motions are set as shown in slice B1–B2. Just beyond the epiglottis region, the flow again accelerates due to reduction in cross-sectional area and a clear high velocity zone develops on the posterior side of sections C1–C2. At the end of the pharynx, a step on the posterior side guides the flow towards the anterior side of the trachea in the form of a laryngeal jet. As a result of this laryngeal jet, two distinct recirculation zones originate at the posterior side at sections E1–E2 and move towards the center as the flow moves further downstream (sections F1–F2). The normalized velocity magnitude contours in Fig. 7 illustrate that the shape of laryngeal jet for LES is much closer to what is observed in experiments, and RANS predicts a more pronounced and longer laryngeal jet compared to LES. It also shows that RANS leads to greater velocities at the anterior tracheal wall, as previously observed at discrete model cross-sections (Fig. 4(b–d)). 5.2. Particle phase Fig. 8 summarizes the simulated total deposition percentages for different particle diameters along with the experimental curve-fit of Grgic et al. (2004), obtained from deposition measurements in seven different model casts representative of over 80 image-based mouth–throat structures. For the 2 and 4 m particles, LES and DES show particle depositions that are much closer to the experimental curve than those obtained with RANS k. (without near-wall anisotropic correction), while for the 8 and 10 m particles, RANS, LES and DES perform equally well.

872

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

Fig. 7. Time-averaged 2-component normalized velocity magnitude contour (ux and uz ) in (a) experiments; (b) LES; (c) RANS.

Alternatively, RANS k. with mean flow tracking, i.e., without EIM, consistently underestimates deposition for all particle diameters greater than 2 m. The same is true for the frozen LES method, possibly because the under-relaxation parameter settings, which yielded good results for LES (Fig. 8), are too dissipative for the frozen LES method. Considering that 5 m is generally referred to as the upper limit of the respirable range for inhalation drugs (represented by the dash-dotted line in Fig. 8), our findings suggest that in the mouth–throat geometry, the prediction of medication aerosol deposition inhaled at normal flow rates were more accurate for LES and DES than for the RANS k. model considered here. At a first glance, DES can then be seen as the preferred method over LES due to reduced computational requirements. To be certain regarding DES being better than LES on 1.2 × 106 mesh points, an additional LES was performed on the DES mesh. It is observed that LES did as good as DES for both fluid and particle phase. This clearly means that LES would remain the preferred method among the models tested, as it obviates solving an additional equation for t required for DES. For the description of particle transport with diameters above 5 m (e.g., in the upper range of air pollutant particle distributions) or for small diameters but inhaled at greater inspiratory flows (e.g., dry powder inhalers), RANS with its vastly lower computational requirements suffices to adequately predict aerosol deposition. In order to better understand the discrepancy in particle deposition between LES (Smagorinsky) and RANS (k.), deposition in the three model sub-parts are shown for the five particle diameters (Fig. 9). It is seen that RANS overestimates larynx/trachea deposition, showing relatively greater discrepancy with LES for the smaller particles. This can be explained by the profound and longer laryngeal jet in case of RANS compared to LES (Fig. 7) and the velocity overestimation near the anterior airway wall (Fig. 4(b–d)). Even though LES and DES are more accurate among the models considered, they pose a limitation when attempting to simulate transient inhalation/exhalation cycles, partly because of the small sampling time requirements and partly due to relatively high time-cost associated with typical transient waveforms. For example, to perform LES/DES simulation of a 2 s transient inhalation waveform as considered by Li et al. (2007) would require ∼ 5 times more computational time compared to the steady inhalation case simulated in the present work. Alternative RANS models may therefore be worth considering. For instance, Matida et al. (2004) applied the near-wall anisotropic corrections to a RANS k. model and showed considerable improvement in deposition over the one obtained with isotropic assumption. However, the fact that the basic fluid flow field is not accurately predicted by RANS (Figs. 4, 5, 7) is a problem. Another possibility is to consider a Reynolds stress model (RSM) which implicitly accounts for near-wall anisotropy, and has been reported

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

873

100 Experimental fit (Grgic 2004) RANS−EIM RANS−MEAN LES Frozen LES DES

90

Deposition Efficiency (%)

80 70 60 50 40 30 20 10 0 10−2

10−1

100

101

Stk.Re0.37 Fig. 8. Simulated total deposition (expressed as % of particles at model inlet) as a function of Stokes and Reynolds number as defined by Grgic et al. (2004). The solid line represents the experimental best fit curve. The dash-dotted line corresponds to a 5 m particle at 30 l/min. In case of RANS, ‘+’ represents turbulent tracking i.e., considering EIM; ‘×’ represents mean flow tracking i.e., without EIM.

70 trachea + larynx pharynx

60

Deposition Efficiency (%)

mouth 50

40

30

20

RANS

10 LES 0

2

4

6 Diameter (µm)

8

10

Fig. 9. Simulated deposition values (expressed as % of particles at model inlet) in three model subparts. Mouth deposition excludes deposition in the inlet tube.

to outperform RANS k. model when predicting nano- and micro-particle deposition in a duct flow (Tian & Ahmadi, 2007). However, as pointed out by these authors, the performance of the different modalities of the RSM in a more complex geometry warrants further investigation.

874

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

6. Conclusions The main conclusions of the present work are summarized as follows: (1) Flow features for the normal breathing condition (30 l/min) were measured in the central sagittal plane of a human mouth–throat model using PIV technique with an aim of providing quality data for the detailed evaluation of CFD simulations in exactly the same mouth–throat geometry. (2) RANS, LES and DES simulations showed that for the fluid phase, LES/DES lead to similar velocity and kineticenergy profiles and all of which compare well with experimental data, as opposed to RANS. (3) Five m is generally referred to as the upper limit of respirable range for inhalation drugs. For the micro-particles below 5 m considered in the present study, LES and DES more closely match experimental aerosol deposition than RANS (without near-wall correction). Simulating LES on the DES mesh showed that LES would be the preferred method among the models tested. (4) For the simulation of particles above 5 m at normal breathing condition, the similar performance of RANS and LES/DES makes RANS the preferred method. Acknowledgments Part of this research was funded by the VUB Research Council in the framework of a horizontal research activity and this funding is gratefully acknowledged. References Breuer, M. (1998). Large eddy simulation of the sub-critical flow past a circular cylinder: Numerical and modeling aspects. International Journal of Numerical Methods in Fluids, 28, 1281–1302. Breuer, M. (2000). A challenging test case for large eddy simulation: High Reynolds number circular cylinder flow. International Journal of Heat and Fluid Flow, 21(5), 648–654. Breuer, M., Baytekin, H. T., & Matida, E. A. (2006). Prediction of aerosol deposition in 90 degree bends using LES and an efficient Lagrangian tracking method. Journal of Aerosol Science, 37, 1407–1428. Breuer, M., Matida, E., & Delgado, A. (2007). Prediction of aerosol drug deposition using Eulerian–Lagrangian method based on LES. In International conference on multiphase flow, July 9–13, Leipzig, Germany. Brouns, M., Jayaraju, S., Lacor, C., Mey, J., Noppen, M., Vincken, W., & Verbanck, S. (2007). Tracheal stenosis: A flow dynamics study. Journal of Applied Physiology, 102, 1178–1184. Corcoran, T. E., & Chigier, N. (2000). Characterization of the laryngeal jet using phase doppler interferometry. Journal of Aerosol Medicine, 13, 125–137. DeHaan, W., & Finlay, W. (2001). In vitro monodisperse aerosol deposition in a mouth and throat with six different inhalation devices. Journal of Aerosol Medicine, 14, 361–367. DeHaan, W., & Finlay, W. (2004). Predicting extrathoracic deposition from dry powder inhalers. Journal of Aerosol Science, 35, 309–331. Elghobashi, S. (1994). On predicting particle-laden turbulent flows. Applied Scientific Research, 52, 309–329. Gemci, T., Corcoran, T., & Chigier, N. (2002). A numerical and experimental study of spray dynamics in a simple mouth throat model. Journal of Aerosol Science, 36, 18–38. Gosman, A. D., & Ioannides, E. (1981). Aspects of computer simulation of liquid-fuelled combustor. AIAA Journal, 81-0323. Grgic, B., Finlay, W. H., Burnell, P. K. P., & Heenan, A. F. (2004). In vitro intersubject and intrasubject deposition measurements in realistic mouth–throat geometries. Journal of Aerosol Science, 35, 1025–1040. Heenan, A. F., Finlay, W. H., Grgic, B., Pollard, A., & Burnell, P. K. P. (2004). An investigation of the relationship between the flow field and regional deposition in realistic extra-thoracic airways. Journal of Aerosol Science, 35, 1013–1023. Heenan, A. F., Matida, E., Pollard, A., & Finlay, W. H. (2003). Experimental measurements and computational modeling of the flow field in an idealized human oropharynx. Experiments in Fluids, 35, 70–84. Hopkins, L., Kelly, J., Wexler, A., & Prasad, A. (2000). Particle image velocimetry measurements in complex geometries. Experiments in Fluids, 29, 91–95. Ilie, M., Matida, E., & Finlay, W. (2008). Asymmetrical aerosol deposition in an idealized mouth with a DPI mouthpiece inlet. Aerosol Science and Technology, 42, 10–17. Jin, H. H., Fan, J. R., Zeng, M. J., & Cen, K. F. (2007). Large eddy simulation of inhaled particle deposition within the human upper respiratory tract. Journal of Aerosol Science, 19, 257–268. Johnstone, A., Heenan, A., Uddin, M., Pollard, M., & Finlay, W. H. (2004). The flow inside a idealized form of the human extra-thoracic airway. Experiments in Fluids, 37, 673–689. Kleinstreuer, C., & Zhang, Z. (2003). Laminar-to-turbulent fluid-particle flows in a human airway model. International Journal of Multiphase Flow, 29, 271–289.

S.T. Jayaraju et al. / Aerosol Science 39 (2008) 862 – 875

875

Li, Z., Kleinstreuer, C., & Zhang, Z. (2007). Simulation of airflow fields and microparticle deposition in realistic human lung airway models. Part I: Airflow patterns. European Journal of Mechanics B/Fluids, 26, 632–649. Luo, X. Y., Hinton, L. S., Liew, T. T., & Tan, K. K. (2004). LES modelling of flow in a simple airway model. Medical Engineering and Physics, 26, 403–413. Matida, E. A., Finlay, W. H., Breuer, M., & Lange, C. F. (2006). Improving prediction of aerosol deposition in an idealized mouth using large eddy simulation. Journal of Aerosol Medicine, 19, 290–300. Matida, E. A., Finlay, W. H., Lange, C. F., & Grgic, B. (2004). Improved numerical simulation of aerosol deposition in an idealized mouth–throat. Journal of Aerosol Science, 35, 1–19. Menter, F. R. (1994). Two equation eddy-viscosity turbulence models for engineering applications. AIAA Journal, 32(8), 1598–1605. Morsi, S., & Alexander, A. (1972). An investigation of particle trajectories in two-phase flow systems. Journal of Fluid Mechanics, 55, 193–208. Scarano, F., & Riethmuller, M. (2000). Advances in iterative multigrid PIV image processing. Experiments in Fluids, 29, 51–60. Stapleton, K. W., Guentsch, E., Hoskinson, M. K., & Finlay, W. H. (2000). On the suitability of k . turbulence modeling for aerosol deposition in the mouth and throat: A comparison with experiment. Journal of Aerosol Science, 31, 739–749. Tian, L., & Ahmadi, G. (2007). Particle deposition in turbulent duct flows—comparison of different model predictions. Journal of Aerosol Science, 38, 377–397. Vinchurkar, S., & Longest, P. (2008). Evaluation of hexahedral, prismatic and hybrid mesh styles for simulating respiratory aerosol dynamics. Computers and Fluids, 37, 317–331. Westerweel, J., & Scarano, F. (2005). Universal outlier detection for PIV data. Experiments in Fluids, 39, 1096–1100. Xi, J., & Longest, P. W. (2007). Transport and deposition of micro-aerosols in realistic and simplified models of the oral airway. Annals of Biomedical Engineering, 35, 560–581.