High resolution turbulence modelling of airflow in an idealised human extra-thoracic airway

Available online at www.sciencedirect.com

Computers & Fluids 37 (2008) 943–964 www.elsevier.com/locate/compfluid

High resolution turbulence modelling of airflow in an idealised human extra-thoracic airway C.G. Ball, M. Uddin, A. Pollard

*

Department of Mechanical and Materials Engineering, Queen’s University, Kingston, Ontario, Canada K7L 3N6 Received 6 October 2006; received in revised form 13 May 2007; accepted 27 July 2007 Available online 6 November 2007

Abstract Computational fluid dynamic (CFD) studies of the flow inside a modelled human extra-thoracic airway (ETA) were conducted to evaluate the performance of several turbulence models in predicting flow inside this complex geometry. Veracity of the computational models is assessed for physiologically accurate flow rates of 10, 15, and 30 l/min by comparison of numerical results with hot-wire [Johnstone A, Uddin M, Pollard A, Heenan A, Finlay WH. The flow inside an idealised form of the human extra-thoracic airway. Exp Fluids 2004;37(5):673–89] and particle image velocimetry (PIV) [Heenan AF, Matida E, Pollard A, Finlay WH. Experimental measurements and computational modelling of the flow field in an idealised extra-thoracic airway. Exp Fluids 2003;35:70–84] mean velocity data for the central sagittal plane of the ETA. Furthermore, flow features predicted by numerical models are compared to those from experimental flow-visualisation studies [Johnstone et al., 2004]. The flow in the ETA is shown to be highly three-dimensional, having strong secondary flows. Ó 2007 Elsevier Ltd. All rights reserved.

1. Introduction Inhaled pharmaceuticals, because of their convenience and rapid action, are widely used for treatment of the most prevalent types of lung disease such as asthma, chronic obstructive pulmonary disease, cystic fibrosis, and others. Anecdotal evidence suggests that both respiratory disease sufferers and their physicians recognise, through negative side-effects, the difficulty in achieving optimal delivery of treatment medicines to the lungs; patients are given specific instructions for best use of their inhalers. Yet despite best medical practise and careful engineering of aerosols, drug delivery rates to the alveolar regions of the lung are minimal; typically 80–95% of an inhaled aerosol dose is deposited in the extra-thoracic airway [3], often with non-benign side-effects, such as inflammation, as the result. Because *

Corresponding author. E-mail addresses: [email protected] (C.G. Ball), uddin@me. queensu.ca (M. Uddin), [email protected] (A. Pollard). URL: http://me.queensu.ca/people/pollard/ (A. Pollard). 0045-7930/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2007.07.021

particle delivery is affected by the flow in the extra-thoracic airway (ETA), it is expected that by modifying the flow field, targetted delivery can be achieved. The ability to predict and control the related phenomena of particle deposition in the ETA (for the anesthetist) and particle delivery to the lungs (for the respirologist) would be a significant improvement over the current delivery patterns of existing inhalers. Although significant effort has been directed toward numerical modelling of deposition in the human airway, accurate prediction of regional drug deposition patterns in the ETA has yet to be achieved. Even for the simplest cases, particle deposition models, such as the Eddy Interaction Model proposed by Hutchinson et al. [4,5], are known to perform poorly. However, inadequacy of particle deposition predictions cannot be solely attributed to the particle tracking models because decoupling of particle dynamics from computational fluid dynamics for verification is incomplete. It is unreasonable to expect accurate numerical predictions of particle deposition to be achieved before the ETA flow field is fully understood, and before single-phase

944

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

CFD models are verified. Nonetheless, attempts to predict particle deposition in the human upper airways are reported in the archival literature. In a comparison study, Zhang and Kleinstreuer [6] undertook to assess the capacity of several RANS models to predict transitional flows in a biological conduit. Four turbulence models were verified against experimental results for an axi-symmetric constricted pipe. For flow predictions in their oral-airway model, which resembles a bent pipe with constrictions (meant to approximate the naso– oro-pharynx, and the larynx and glottis), the results from four simulations using different numerical models were compared to each other. The validity of the numerical predictions of flow in the simple upper airway was not assessed by comparison with experimental results. The analysis presented in [6,7] appears to be a precursor to particle deposition studies for the same oral-airway model as in [7–9] and others. Regardless, there is very little similarity between the simple airway of Zhang et al. and the idealised ETA used in the work described herein. In other numerical deposition studies, such as [10,11], and others by the same group, particle deposition rates (as percentages of total administered dose) rather than deposition patterns are presented. Furthermore, validation of the CFD airflow solution is limited to a single low-resolution study of the performance of the k– model. To accurately predict the flow in the complex geometry of the ETA is a daunting challenge. Comprising the oral and nasal cavities, the pharynx, the larynx, and the trachea (as shown in Fig. 1), the ETA presents a variety of crosssectional shapes and diameters (as shown in Fig. 2 for

Fig. 2. The idealised ETA, feature names, cross-sectional shapes, and model dimensions (in mm).

the idealised ETA) so that the flow can range from laminar-like to turbulent. For the flow rates considered, the local Reynolds number, based on local mean diameter and local bulk velocity (defined in Section 2.3), ranges from approximately 480 to over 3100 (as shown in Table 1). By comparison with the PIV data of Heenan et al. [2] and with the hot-wire measurement and flow-visualisation studies of Johnstone et al. [1], the veracity of several Reynolds Averaged Navier–Stokes (RANS) turbulence models to predict the complex flow of the ETA is assessed in this paper. While the accuracy of numerical predictions is assessed by comparison with experimental results, it is recognised that, unfortunately, there exists some question as to the best method (experimental or numerical) for verification of results. Given the complex geometry, which approach is the best arbiter? The work of Heenan et al. [2] and Johnstone et al. [1] demonstrate that complete trust in PIV and hot-wire measurements cannot be made – hot-wire measurements are unreliable in regions with recirculation and secondary flow phenomena, PIV measurements are compromised by optical access and seeding limitations. Even though the strengths and limitations of the various turbulence models are well known [12–15], they nonetheless cast doubt onto the veracity of predicted results. Traditional Table 1 Reynolds number and bulk velocity by station Flow rate (l/min)

Fig. 1. The human airways [37].

10

15

Section

UB

Re

UB

Inlet 1 2 3 4 5 6 7

0.587 0.291 0.275 0.542 1.178 0.623 1.367 0.778

677 477 464 651 960 698 1034 780

0.880 0.436 0.412 0.813 1.767 0.934 2.050 1.167

30 Re

UB

1016 1.760 715 0.873 695 0.825 976 1.625 1439 3.534 1047 1.870 1551 4.101 1170 2.334 pffiffiffi Bulk velocity UB in m/s, Reynolds number Re  U B A=m.

Re 2032 1431 1391 1953 2879 2094 3101 2340

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

direct numerical simulation (DNS) techniques are not suited to handle the complicated geometry of the ETA. With knowledge of the characteristics of these multiple flow elucidation methods, the best from each are taken to build a complete understanding of the ETA flow. Thorough and accurate understanding of the complicated ETA flows can be achieved by concurrent evaluation of experimental and computational results. This work is an evaluation of the performance of RANS turbulence models applied to an idealised human ETA, and a presentation of the types of flow patterns developed as air passes from the mouth to the trachea. 2. Methods 2.1. The idealised ETA geometry Inter-subject variability in the human ETA, as can be seen in Heenan et al. [16] and documented in [17,18], presents a significant obstacle to the generation of meaningful findings that can be applied to the general population. The physical model chosen for this work, and for the experimental investigations of Heenan et al. [2] and Johnstone et al. [1], is a physiologically accurate ‘‘average’’ geometry of the human upper airway. It is a variation of that described in [2,18], whose geometry is based on data from computed tomography (CT) scans, magnetic resonance imaging (MRI) scans, direct observation of human airways, and data from archival literature. As can be seen from comparison of the model with an actual ETA, as shown in Figs. 1 and 3, it comprises all the fundamental features (such as the epiglottis, larynx, etc.) of the ETA. The complexity of the geometry is revealed by examination of its cross-sections, as can be seen in Fig. 2. 2.2. Meshing, grid resolution, and turbulence models A CAD version of the idealised ETA generated in Pro/ ENGINEER (PTC, Needham, MA) was exported in Initial

Graphics Exchange Specification (IGES) format for meshing in GAMBIT (v.2.0.4, Fluent Inc., Lebanon, NH). The tetrahedral-hybrid meshing scheme was applied uniformly throughout the geometry. Although wall-biased mesh refinements are customarily employed in wall-bounded flow simulations, due to limited a priori knowledge of the ETA flow patterns, no such meshing restrictions were imposed in this study. Experimental evidence suggests that deficient regional refinement of the mesh could result in failure to capture significant features of the flow because it is not yet fully understood where such refinements should be applied (in addition to near the walls). The grid size, or mesh refinement, varies with inhalation flow rate to maintain a cell Reynolds number in the range 5 < Recell < 85. Cell Reynolds number is defined as Recell  Q/Lcellm, where Lcell is the control volume size. The grid for the 30 l/min case had over nine million control volumes, with a nominal cell size of 0.0325 mm. This grid density is 40 times that used by Stapleton et al. [18] and 10 times that used by Heenan et al. [2]. The wall-bound control volume size is y+ = O(1) (this is important for the k–x simulations where wall functions are not used). The maximum y+  2.8, found at the posterior wall near station 4 in the 30 l/min case, is still sufficiently small to capture wall effects. Turbulence was modelled by implementing various modules of the commercial CFD code FLUENT (v.6.1, Fluent Inc., Lebanon, NH). The more commonly used turbulence models were compared in this study; specifically, these are the Standard and the Shear Stress Transport (SST) versions of the k–x model, the Standard and the Renormalisation Group (RNG) versions of the k– model, and the Reynolds stress model (RSM). Whenever possible, wall functions were not implemented (that is, for the k–x simulations); a more accurate solution is expected by integrating to the wall rather than making an assumption (i.e., the wall function) regarding the near wall behaviour of the flow [12]. The standard wall function was used for the k– and RSM turbulence calculations. For the all models, pressure was solved by a first-order discretisation scheme in which pressure values at control volume faces are interpolated using momentum equation coefficients [19]. Pressure–velocity coupling was accomplished using the SIMPLE scheme. All other flow variables (momentum, k, , x, and the Reynolds stresses) were discretised using the QUICK scheme. A brief description of the various models of turbulence is in order. The time-averaged continuity and momentum equations for three-dimensional incompressible steady state flows are, respectively, oðqui Þ ¼0 oxi     oðqui uj Þ op o oui ouj 0 0 ¼ þ l þ  qui uj oxj oxi oxj oxj oxi

Fig. 3. Sagittal section of ETA model, showing areas of expected significant flow phenomena (a–i).

945

ð1Þ ð2Þ

From the Reynolds decomposition, the second-order correlations, Rij ¼ qu0i u0j , that result are known as the Reynolds stresses.

946

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

2.2.1. The standard k– model For the standard k– model, turbulent kinetic energy k and its dissipation rate  are related by mt = Clk2/. The transport equations for k and  are given by    o o lt ok ðqui kÞ ¼ lþ ð3Þ þ Gk  q þ S k oxi oxi rk oxi    o o l o  2 ðqui Þ ¼ lþ t þ C 1 Gk  C 2 q þ S  ð4Þ oxi oxi k r oxi k In these equations, Gk is the generation of turbulence kinetic energy due to mean velocity gradients; C1, C2, and Cl are constants rk and r are the Prandtl numbers for k and , respectively; Sk and S are source terms. The values of the constants used were determined from fundamental shear flows in air and water [20,21] C u ¼ 0:09; r ¼ 1:3

C 1 ¼ 1:44;

C 2 ¼ 1:92;

rk ¼ 1:0;

2.2.2. The renormalisation group k– model Mathematical ‘‘renormalisation group’’ (RNG) methods were used to derive the RNG k– turbulence model from the instantaneous Navier–Stokes equations. The result is a model with different constants, and additional terms and functions in the transport equations for k and . The transport equations for the RNG k– model are written:   o o ok ðqui kÞ ¼ ak leff ð5Þ þ Gk  q þ S k oxi oxi oxj   o o ok  2 ðqui Þ ¼ ak leff þ C 1 G k  C 2 q  R þ S  oxi oxi oxj k k ð6Þ For the RNG k– model: turbulence production Gk is modelled in the same manner as for the standard k– model; ak and a are the inverse effective Prandtl numbers for k and , respectively; the model constants C1, C2, and Cl are derived using RNG theory. The differential equation for turbulent viscosity in the RNG method is given by  2  ^m qk d pffiffiffiffiffi ¼ 1:72 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi d^m ð7Þ 3 l ^m  1 þ C m where ^m ¼ leff =l and Cm  100. The additional R term is given by

2.2.3. The k–x model The model as proposed by Wilcox [22]; its governing equations are   o o ok ðqui kÞ ¼ Ck ð9Þ þ qG  Y k þ S k oxi oxi oxi   o o ox ðqui kÞ ¼ Cx ð10Þ þ Gx  Y x þ S x oxi oxi oxi where the turbulent viscosity is computed: lt ¼ a q

k x

ð11Þ

In the k–x transport equations, Gk represents generation of turbulence kinetic energy by mean velocity gradients; Gx represents generation of x; the terms Ck and Cx are the effective diffusivities of k and x, respectively; Yk and Yx are the dissipation of k and x, respectively; Sk and Sx are source terms. The effective diffusivities are Ck = l + lt/rk and Cx = l + lt/rx where rk and rx are the turbulent Prandtl numbers. A low-Reynolds number correction for turbulent viscosity is included with the coefficient a*, given by    a0 þ Ret =Rk   a ¼ a1 ð12Þ 1 þ Ret =Rk where Ret = qk/lx and Rk ¼ 6;

Rx ¼ 2:95;

a0 ¼

bi ; 3

bi ¼ 0:072

In the high-Reynolds form a ¼ a1 ¼ 1. The constants used for this implementation of the k–x model are rk ¼ 2;

rx ¼ 2

b 1 a0 ¼ ; a0 ¼ 3 9 b ¼ 0:072; Reb ¼ 8 Rek ¼ 6;

Rex ¼ 2:95

And, the turbulent Reynolds number, Ret is calculated using Ret ¼

qk xl

ð13Þ

The remaining constants C1 and C2 were derived analytically to be

2.2.4. The shear stress transport k–x model In the SST k–x model, the formulation for turbulent viscosity is modified to include turbulent shear stress transport. Further, a cross-diffusion term in the x equation and a blending function are incorporated to improve model performance in both the near-wall and in the far-field regions. The transport equations take the form   o o ok ðqkui Þ ¼ Ck ð14Þ þ Gk  Y k þ S k oxi oxj oxj

C 1 ¼ 1:42;

and

3

R ¼

2

C l qm ð1  m=m0 Þ  1 þ bm3 k

ð8Þ

where m

Sk ; 

m0 ¼ 4:38;

b ¼ 0:012

C 2 ¼ 1:68;

C l ¼ 0:0845

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

  o o ox ðqxui Þ ¼ Cx þ Gx  Y x þ Dx þ S x oxi oxj oxj

ð15Þ

Here, Gk, Yk, Yx representing turbulence kinetic energy, x generation, dissipation of k, dissipation of x, respectively, are all calculated in the same manner as for the standard k–x model. Calculation of turbulence production of x, Gx, the effective diffusivities, Ck and Cx, and the cross-diffusion term, Dx are accomplished somewhat differently. Turbulent viscosity is computed by qk 1 ð16Þ  x maxð1=a ; XF 2 =a1 xÞ pffiffiffiffiffiffiffiffiffiffiffiffiffiffi where X  2Xij Xij (Xij is the mean rate-of-rotation tensor, a* is as in Eq. (12)). F2 is a blending function. Except for those listed here, the constants for the SST and for the standard k–x models take the same values: a1 ¼ 0:31; bi;1 ¼ 0:075; bi;2 ¼ 0:0828

lt ¼

rk;1 ¼ 1:176;

rk;2 ¼ 1:0;

rx;1 ¼ 2:0;

rx;2 ¼ 1:168

2.2.5. The Reynolds stress model In the Reynolds stress model, the individual Reynolds stresses are computed directly using differential transport equations and used to close the Reynolds-averaged momentum equation. The transport equation of the Reynolds stresses is written as o o ðqu0i u0j Þ þ ðquk u0i u0j Þ ot ox k |fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} Local Time Derivative

C ij Convection

  o o o 0 0 0 0 0 0 0 ¼ ½qu1 uj uk þ pðdkj ui þ dik uj Þ þ l ðu u Þ oxk oxk oxk i j |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} DT ;ij Turbulent Diffusion

947

imation of the inlet condition as steady is justified, as noted by Chang and Masry [23] who show that, based on quasisteady flow criteria, the physiological flow regime that occurs during the major part of the inhalation cycle is well matched by this assumption. As a result of flow visualisation studies [1,24] and ongoing numerical work, it is recognised that in some regions of the ETA the flow may be unsteady. The viscous forces of a flow are related pffiffiffiffiffiffiffiffiffiffiffito ffi the unsteadiness by the Womersley number, a ¼ Re=St (St is the Strouhal number). In general, the unsteadiness can result from either of two sources: the unsteady inhalation (similar to the original sinusoidal inlet condition used by Womersley in his solution to oscillatory laminar flow in a pipe [25]), and/or secondary flow unsteadiness within the ETA [26]. In this work, steady inflow conditions were used; from this perspective, the Womersley number is inapplicable. However, there may be a time scale associated with the internal flow dynamics (such as vortex shedding) that could be related to a Strouhal number, St. In keeping with the foundation works of Heenan et al. [2] and Johnstone et al. [1], the numerical simulations of the current study also employed the assumption that the inspiratory ETA flow is time-steady. A uniform velocity profile was imposed at the inlet plane of a pipe extension to the oral cavity. In the experiments, a bell-mouth contraction was used to produce a uniform inlet profile at the entrance to the oral cavity; however, tests by Heenan et al. [2] show no effect of varying inlet condition. An outflow boundary condition was specified at the outlet. The walls of the model were specified as being dry and inflexible. This seemingly significant modelling assumption is justified because (1) at physiological flow rates, there is

DL;ij Molecular Diffusion

  ouj oui 0 0 0 0  q ui uk þ ui uk  qbðgi u0j h þ gj u0i hÞ oxk oxk |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Gij Buoyancy Production 

P ij Stress Production

 ou0i ou0j ou0 ou0j þp þ  2l i oxj oxi oxk oxk |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflffl{zfflfflfflfflfflffl} /ij Pressure–Strain

ij Dissipation

 2qXk ðu0j u0m ikm þ u0i u0m jkm Þ þ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} F ij Production by System Rotation

S |{z}

ð17Þ

Source Term

These are then modelled. The constants used for the current implementation of the RSM are C l ¼ 0:09; C 1 ¼ 1:44; C 2 ¼ 1:92 C 1ps ¼ 1:8;

C 2ps ¼ 0:6;

C 01ps ¼ 0:5;

C 02ps ¼ 0:3

Turbulent viscosity is computed using lt = qCll2/. The C terms are used in computing turbulence dissipation. The Cps terms are used in modelling pressure–strain. 2.3. Boundary conditions and simulation parameters The simulations were carried out for physiologically accurate steady flow rates of 10, 15, and 30 l/min. Approx-

Fig. 4. ETA inlet at station 1, showing alignment of local co-ordinate axes.

948

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

no significant difference in the viscous pressure drop caused by a mucous lining as compared to that caused by a simple constriction [27] and (2) length and diameter variations of the upper airways are small, the change rates of these variations are negligible as compared to airflow velocities [28]. Finally, the working fluid is modelled as atmospheric air. Because of the varying cross-sectional geometry of the ETA, the pffiffiffi Reynolds number was calculated as Re  Q= Am, where Q is the volumetric flow rate and A

is the local cross-sectional area, with density assumed to be constant. The bulk velocity was defined as Ubulk  Q/A. For each flow rate, the bulk velocity and Reynolds number for the inlet, outlet, and each measurement station are given in Table 1. Transitional flows are indicated by these ranges of Reynolds number; it was confirmed that the code produces the same parabolic velocity profile when either a k–x or a laminar pipe flow model was considered (Re = 1000). In order to be consistent with previous ETA

3

2.5

2

K ⎯ 1.5 UB

JUPHF Single-Wire U (2004) ⎯⎯⎯ JUPHF X-Wire √U +V (2004) Standard k - ω SST k - ω Standard k - ε RNG k - ε RSM

1

0.5

0

0

0.5

y ⎯ ⎯√Α

1

Fig. 5. Station 1: 30 l/min velocity magnitude profiles.

3.0

2.5

2.0

1.5

Ux ⎯ UB

Standard k - ε Standard k - ω RNG k - ε SST k - ω RSM

1.0

0.5

0.0 0 -0.5

0.5

y ⎯ ⎯√Α

1

Fig. 6. Station 1: Mean streamwise velocity profiles normalised by bulk velocity, 30 l/min.

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

studies and ongoing projects, the results are rescaled and presented in terms of the 1:1 ETA prototype. 2.4. Co-ordinate systems The use of a single global reference frame would, because of the irregular geometry and station placement, make interpretation of flow features difficult and ambiguous. Local co-ordinate systems, corresponding to each of the measurement stations were therefore defined. Demonstrated using station 1 in Fig. 4, the y-axis is along the wall-normal station traverse, with y = 0 being on the superior (for stations 1–3) or posterior (for stations 4–7) wall of the model; the positive x-axis is aligned with the bulk flow direction of the station; and following the right-hand rule, z is perpendicular to the sagittal plane of the geometry. Uppercase U, V, and W denote mean velocity components in the x-, y-, and z-directions, respectively. Overall velocity pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi magnitude is given by K ¼ U 2 þ V 2 þ Wp2 .ffiffiffiLengths and velocities were non-dimensionalised by A and Ubulk, respectively.

949

2.5. Computing Simulations were run at the High Performance Computing Virtual Laboratory (HPCVL) at Queen’s University using Sunfire 1.05 GHz UltraSparc III Cu processors. Performance tests indicated maximum scalability to eight processors; cases were partitioned and parallelised to four CPUs using Fluent’s autopartitioning scheme. Wall-clock simulation time for the 30 l/min case was typically 2–3 days for the 2-equation models, and 10–12 days for the Reynolds stress model. Calculations ran until the model variables reached a convergence level of O(104); continuity, momentum, and the turbulence quantities (k, , x, and the stresses) were monitored. 3. Results and discussion The experimental data were obtained in our laboratory using a standard technique called hot-wire anemometry.

Fig. 8. Station 2: Direction of probe traverse.

Fig. 7. Oral cavity: Velocity vectors coloured by velocity magnitude, standard k–x 10 l/min. Top: showing the full range of velocities in the oral cavity, from 0 to 1.1 m/s; bottom: the colour spectrum is limited to the velocity range 0–0.2 m/s.

Fig. 9. Station 2, frontal plane: Velocity vectors coloured by velocity magnitude, RSM 10 l/min.

950

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

This method uses very fine wires (5 lm diameter, 200:1 length to diameter ratio). Thus, while data are stated to be obtained at a particular traverse location, in fact, each measurement value is spatially averaged. Also, a single hot-wire can neither distinguish flow direction nor velocity components. Indeed, the time-averaged output represents the absolute value of the local flow vector, jKj.

3.1. Station 1 Station 1 is located immediately downstream of the inlet to the oral cavity, as shown in Fig. 4. The upper curvature of the mouth causes flow deceleration, an increase in pressure, and the formation of a recirculation bubble. The lower region, comprising the backward facing step (BFS)

3 JUPHF Single-Wire U (2004) ⎯⎯⎯ JUPHF X-Wire √U +V (2004) Standard k - ω SST k - ω Standard k - ε RNG k - ε RSM

2.5

2

K ⎯ 1.5 UB

1

0.5

0

0

0.2

0.4

0.6

0.8

1

y ⎯ ⎯√Α Fig. 10. Station 2: 15 l/min velocity magnitude profiles.

2.5

JUPHF Single-Wire U (2004) 5mm upstream Station 2 5 mm downstream

2.0

1.5

K ⎯ UB 1.0

0.5

0.0

0

0.2

0.4

y ⎯ ⎯√Α

0.6

0.8

Fig. 11. Station 2: Mean velocity magnitude profiles normalised by bulk velocity; k–x, 15 l/min. Effect of translation of the measurement station.

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

of the lower teeth, causes separation and recirculation. Shortly, downstream of station 1, the oral cavity expands laterally; this creates instabilities in the first regions of the ETA during inhalation that may be carried downstream to grow or to be damped out, depending on Reynolds number. Being near to the round-pipe inlet, the flow features at this station are similar to those of flow in a circular pipe; and so, all of the turbulence models are expected to pro-

951

duce acceptable flow predictions here. Computational and experimental velocity magnitude results for the 30 l/ min case are shown in Fig. 5. As compared to experiment, the k– models, with incorrect sloped ‘‘tops’’ and poor match near x  0.3, produced inferior results to those predicted by the k–x and Reynolds stress models. It is expected that under these conditions the k– model will exhibit poor performance [12]. There is some difficulty with the predictions using the Reynolds stress model,

2.5

JUPHF Single-Wire U (2004) Station 2: Standard k - ω Station 2 +10 degrees Station 2 +5 degrees Station 2 -5 degrees Station 2 -10 degrees

2.0

1.5

K ⎯ UB 1.0

0.5

0.0

0

0.2

0.4

0.6

0.8

y ⎯ ⎯√Α 2.5

JUPHF Single-Wire U (2004) Station 2: Standard k - ω Station 2 +10 degrees Station 2 +5 degrees Station 2 -5 degrees Station 2 -10 degrees

2.0

1.5

K ⎯ UB 1.0

0.5

0.0

0

0.2

0.4

y ⎯ ⎯√Α

0.6

0.8

Fig. 12. Station 2: Mean velocity magnitude profiles normalised by bulk velocity; k–x, 15 l/min. Top: effect of station-line rotation in the xy-plane; bottom: effect of station-line rotation in the xz-plane.

952

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

particularly for higher flow rates, as first seen in the region pffiffiffi y= A P 1:0 in Fig. 5. In the streamwise velocity plot of Fig. 6, no backflow is predicted by any model in this region; the RSM presents non-physical flow in the y- and z-directions. Downstream predictions by the RSM at higher flow rates are also clearly non-physical; they violate continuity laws and so are not reproduced here.

Evidence of the recirculation bubbles generated near station 1 can be seen pffiffiffi in velocity profiles of Figs. 5 and 6 in the region 0 6 y= A 6 0:15. Similar phenomena are observed at station 5, downstream of another backward facing step; and at station 7, downstream of a rapid expansion. The expected recirculation bubble behind the BFS of the lower teeth is barely discernible in both the velocity profile plots and the velocity magnitude contour plots of the central sagittal plane of the oral cavity, as shown in Fig. 7. This seemingly missing feature appears to be shifted far downstream of station 1; according to the k–x model 30 l/min prediction, it may be found just upstream of station 2. Whether this predicted behaviour is an artefact of the laterally curving BFS of the lower teeth, or whether it is an artefact of errors with the turbulence model is unknown. It is possible that, in addition to the curvature, the laterally varying step height, in conjunction with the laterally expanding oral cavity, produces recirculation phenomena that are minimally, if at all, observed in the central sagittal plane. This is the first of several morphological features which complicate the flow in the ETA. 3.2. Station 2

Fig. 13. Measurement stations 3 and 4, showing direction of probe traverse.

Station 2, as shown in Fig. 8, is also in the oral cavity; however, there is question as to whether it is entirely downstream of the recirculation zones generated by the backward facing step and the upper expansion near the inlet to the oral cavity. The time varying nature of these zones was revealed in flow visualisation studies [1,24]. The unsteady nature of the recirculation bubbles, coupled with secondary flow phenomena generated by the expansion of

2

1.5

K ⎯ UB 1 JUPHF Single-Wire U (2004) ⎯⎯⎯ JUPHF X-Wire √U +V (2004) Standard k - ω SST k - ω Standard k - ε RNG k - ε RSM

0.5

0

0

0.2

0.4

0.6

0.8

y ⎯ ⎯√Α Fig. 14. Station 3: Mean velocity magnitude profiles normalised by bulk velocity, 10 l/min.

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

the oral cavity (as described in [1]) may account for the weak coincidence here between experimental data and results produced by steady-state calculations. Both experimental and numerical results are based on a steady-state assumption, yet they are not entirely alike. The experimental results are compiled from ensemble averaged 5-s measurement sets while the numerical simulations present an all-time steady-state solution; it is this difference that may result in some discrepancy. The secondary vortices are created by flow driven towards the sides and down into the lower-lateral corners of the oral cavity; these phenomena are predicted by the k–x and RSM models used in this study (as shown for the RSM in Fig. 9). The rapid expansion created at the inlet to the oral cavity is the likely source of these vortices [1]. Despite the variability in results, each of the k–x and the Reynolds stress models calculated velocity profiles have a multi-stepped shape reminiscent of the experimental profiles, as shown in Fig. 10. Evidentially, the k– model predictions fail by a large margin to capture the flow physics. The exact nature and three-dimensional extent of the secondary vortices and of recirculation on the floor of the oral cavity is unknown. However, it is a possible contributor to the large discrepancy seen in the results at this station. Velocity magnitude vectors in the central sagittal plane of the oral cavity, as shown in Fig. 7, reveal unusual flow development in the lower half of the posterior region of the oral cavity. This could be incipient separation of the flow, which is a significant challenge to CFD methods [29,30].

953

Fig. 16. The oro-pharynx and pharynx: stations 4 and 5.

The possibility exists that imperfect replication of the experimental measurement station location and the attendant spatial averaging effects are the causes of discrepancy in station 2 results. Both translational adjustment (the computational station is moved 5 mm upstream and downstream, along the local x-axis, see Fig. 11) and angular adjustment (the station is rotated about the local co-ordinate system, about the x-axis and about the z-axis, see Fig. 12) are seen to produce considerable variation in the velocity magnitude profiles. It is surprising to see that despite poor congruence here, computational results compare quite well at downstream measurement stations; this suggests a possible unknown experimental error at station 2.

2

1.5

K ⎯ UB 1

JUPHF Single-Wire U (2004) ⎯⎯⎯ JUPHF X-Wire √U +V (2004) Standard k - ω SST k - ω Standard k - ε RNG k - ε RSM

0.5

0

0

0.1

0.2

0.3

y ⎯ ⎯√Α

0.4

0.5

Fig. 15. Station 4: Mean velocity magnitude profiles normalised by bulk velocity, 10 l/min.

954

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

3.3. Stations 3 and 4 Stations 3 and 4 are located in the oro-pharyngeal region of the idealised ETA, as shown in Fig. 13. At station 3, the ETA passage is characterised by decreasing crosssectional area and by the local curvature; these cause flow acceleration and momentum redistribution. Station 4 is in the region of smallest cross-sectional are of the ETA.

At station 3, the general trends observed experimentally are not predicted by any of the turbulence models, as shown in Fig. 14. The double peak in the mean velocity profile suggests flow acceleration along the superior surface and an emerging wall-jet [31] along the posterior surface. In the ETA, the highest bulk velocity is found at station 4. Here, for the 30 l/min case, Re  3150 and y+ of the posterior wall-bound cells is 2.8, small enough to capture

3 JUPHF Single-Wire U (2004) ⎯⎯⎯ JUPHF X-Wire √U +V (2004) Standard k - ω SST k - ω Standard k - ε RNG k - ε RSM

2.5

2

K ⎯ UB 1.5

1

0.5

0

0

0.2

0.4

0.6

0.8

y ⎯ ⎯√Α Fig. 17. Station 5. Mean velocity magnitude profile normalised by bulk velocity, 10 l/min.

3 JUPHF Single-Wire U (2004) ⎯⎯⎯ JUPHF X-Wire √U +V (2004) Standard k - ω SST k - ω Standard k - ε RNG k - ε RSM

2.5

2

K ⎯ UB 1.5

1

0.5

0

0

0.2

0.4

y ⎯ ⎯√Α

0.6

0.8

Fig. 18. Station 5. Mean velocity magnitude profile normalised by bulk velocity, 15 l/min.

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

wall effects; both of these values are the largest found within the model. As a result of the flow constriction leading up to this station, velocity profiles (as seen both experimentally and computationally) develop a close-to top-hat shape, as seen in Fig. 15. However, both single- and X-wire experimental results show a slight depression in this top-hat shape, a remnant of the double-peaked profile of station 3, near pffiffiffi y= A  0:2 that is not apparent in the numerical results. The apparent premature velocity drop-off of experimental pffiffiffi results near y= A  0:45 at station 4 is believed to be

955

due to a combination of wall effects and placement limitations of the hot-wire probe. Despite the incredible incongruence at station 2, good agreement of experimental and computational results is exhibited at downstream stations. The curvature and constriction near stations 3 and 4 may have the effect of increasing the Reynolds number for transition to turbulence – the secondary flows observed here could explain the transition number increase from Re  2000 for straight pipes to Re  6000 for curved pipes [32]. The implication is that in the oro-pharynx of the ETA, secondary flows

3 JUPHF Single-Wire U (2004) ⎯⎯⎯ JUPHF X-Wire √U +V (2004) Standard k - ω SST k - ω Standard k - ε RNG k - ε

2.5

2

K ⎯ UB 1.5

1

0.5

0

0

0.2

0.4

y ⎯ ⎯√Α

0.6

0.8

Fig. 19. Station 5. Mean velocity magnitude profile normalised by bulk velocity, 30 l/min.

3.0

2.5

2.0

Standard k - ω SST k - ω Standard k - ε RNG k - ε RSM

1.5

Ux ⎯ UB

1.0

0.5

0.0 0.2 -0.5

0.4

0.6

y ⎯ ⎯√Α

Fig. 20. Station 5. Mean velocity magnitude profile normalised by bulk velocity, 15 l/min.

956

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

contribute to the decay of turbulence observed as the flow moves from the upper to the lower airways. It is perhaps because of this, and the restriction leading up to station 4 that the flow is, in fact, reset, and that the turbulence models generate plausible results at downstream stations despite the discrepancy observed at station 2. 3.4. Station 5 and the epiglottis The second backward facing step of the idealised ETA is located at the upstream end of the modelled pharynx, as shown in Fig. 16; it represents the area where the nasoand oro-pharynges merge. Modelling of the naso-pharynx as a step (rather than as an orifice through which air might flow) is deemed acceptable because of the negligible airflow through the nose during mouth breathing [3]. At the step, the cross-sectional area expands by more than 100% from a quasi-half- to a quasi-full-elliptical shape; thus, a jet-like flow (see for example [33,34]) is expected downstream of the larynx. Suffice it to say, to the authors’ knowledge, there exists no systematic study of the flow in this type of complex geometry. Measurement station 5 is downstream of this backward facing step. It intersects the recirculation zone created by the step and the jet-like flow from the constriction at station 4. Experimental and computational results are compared in Figs. 17–19 which show the velocity magnitude profiles for various flow rates. There is moderate discrepancy between computational and experimental results in pffiffiffi the region 0 6 y= A 6 0:15. However, this is a region with significant reverse flow, as can be seen in Figs. 20–22; stationary hot-wire probes are incapable of discriminating such flows. Furthermore, Johnstone et al. [1] observed that at station 5 the velocity vectors fall outside the probe calibration range for thispregion of the cross-section. The hotffiffiffi wire data for 0 6 y= A 6 0:3 at station 5 must thus be

Fig. 21. Pharynx, central sagittal plane: velocity vectors from PIV measurements, taken from [2].

Fig. 22. Pharynx, sagittal plane. Velocity vectors coloured by magnitude, k–x, 10 l/min.

regarded as untrustworthy. The RANS predicted velocity results instead can be compared to PIV results obtained by Heenan et al. [2]. The recirculation in this region can be seen in both the measured and the calculated velocity vectors plots of Figs. 21 and 22, respectively. The PIV measured velocity field is shown in Fig. 21; the k–x turbulence model results are shown in Fig. 22 (the BFS is at the upper right-hand corner). Additionally, strong secondary flows damped by the constriction leading up to station 4 are renewed here. Because of the curvature of the walls, the height of the backward facing step is not constant across its length. These two morphological features are the source of the secondary flows, shown in Fig. 23, that are carried around the epiglottis and into the larynx (as shown in Figs. 22 and 24, and the flow visualisations of [1]). The epiglottis is located downstream of station 5. On the anterior side of the epiglottis, a recirculating flow region is observed. This time and space varying bubble will modify the flow as it passes over the epiglottis and into the larynx. As the flow passes over the leading edge of the epiglottis (which is not planar, like in a standard forward facing step), a shear layer develops which evolves downstream through an almost sinusoidal path to the larynx. Thus, the three-dimensional separation is followed by additional three-dimensional strain imposed by navigating the epiglottal–laryngeal path. This is a significant challenge to both

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

experimental and numerical measurement techniques. A vector plot of the lower pharynx and the epiglottal region, shown in Fig. 22 shows some of the complexity of this region of the flow; it is compared to the flow visualisation image in Fig. 24, taken from [1]. The authors expect that neither the classic backward facing step of the pharynx nor the sharp leading edge modelled epiglottis are true to the smooth and deformable walls of an actual living upper airway. However, recall that implementation of such boundary conditions in a CFD code or in experiment is presently not practicable. Rather, such limitations are accepted; this work is considered a step toward more accurate modelling of physiological systems. 3.5. Stations 6 and 7 The larynx, the location for measurement station 6, is most likely the significant contributor to airflow resistance in the ETA [35]. At this station, shown in Fig. 25, effective-

957

ness of the turbulence models to accurately predict the velocity profiles is similar to that seen at upstream stations. The flow continues to be highly three-dimensional; secondary flows, as generated at the entrance to the pharynx, are modified as a result of flow around and acceleration over the three-dimensional epiglottis. The laryngeal jet, which appears to form due to acceleration around the structures of the larynx, is seen in the upper portion of Fig. 26; additional instabilities in this three-dimensional flow field are seen in the anterior portion of this figure as two pairs of counter-rotating vortices. The laryngeal jet appears to form primarily due to required acceleration around the structures of the larynx, rather than acceleration due to constriction in the larynx. Station 7 is located near the caudal end of the trachea, just downstream of the larynx (as shown in Fig. 25). This station is in a region where the laryngeal jet impinges on the anterior wall of the trachea; this impingement is believed to be a significant contributor to wall impaction

Fig. 23. Station 5, coronal plane. Velocity vectors projected onto the yz-plane, coloured by magnitude, k–x, 10 l/min.

Fig. 24. Pharynx, central sagittal plane: flow visualisation snapshots, taken from [1].

958

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

Fig. 27. Station 7, coronal plane. Velocity vectors coloured by magnitude, RSM, 10 l/min. Fig. 25. Larynx and trachea: Measurement stations 6 and 7. X

Z

Y

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Fig. 28. Station 7, coronal plane. Velocity vectors coloured by magnitude, RSM 10 l/min. Fig. 26. Station 6, coronal plane. Velocity vectors coloured by magnitude, k–x, 10 l/min.

of inhaled pharmaceutical aerosols [2]. This impingement is seen in Figs. 27 and 28 as the red-coloured1 arrows; this is the laryngeal jet as it enters the trachea. 1 For the interpretation of colour in these figures, the reader is referred to the Web version of this article.

The effect of the laryngeal jet is clearly seen in the velocity profiles of measurement station 7, as shown in Figs. 29 and 30. Because of the expansion, the velocity profiles seen here are similar to those of station 5, having a significant recirculation zone in the posterior region near the inlet to the trachea. Also, as at station 5, the hot-wire results are considered pffiffiffi questionable in this recirculation zone, for 0 6 y= A 6 0:6, as shown in Figs. 29 and 30.

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

959

4

3.5

JUPHF Single-Wire U (2004) ⎯⎯⎯ JUPHF X-Wire √U +V (2004) Standard k - ω SST k - ω Standard k - ε RNG k - ε RSM

3

2.5

K ⎯ UB

2

1.5

1

0.5

0

0

0.5

y ⎯ ⎯√Α

1

Fig. 29. Station 7: Mean velocity magnitude profiles normalised by bulk velocity, 10 l/min.

3.0

2.5

2.0

Standard k - ω SST k - ω Standard k - ε RSM

1.5

Ux ⎯ UB 1.0

0.5

0.0 0.5 -0.5

y ⎯ ⎯√Α

1

Fig. 30. Station 7: Mean streamwise velocity normalised by bulk velocity, 10 l/min.

Experimental flow visualisations in this region of the ETA reveal highly swirling turbulent flow [24]; non-distinct streamlines confirm extensive turbulent mixing of the flow. This is also seen in velocity vector plots; however, in this region, as elsewhere in the ETA, the nature of these plots is highly dependant on the turbulence model used. Evidence of this variability inpmodel results is seen in the recirffiffiffi culation region for 0 6 y= A 6 0:6, in Figs. 29 and 30. The

k–x model provides the best agreement with experimental results, as seen by examination of the slopes of the curves near x  1.2 and 0.5 6 x 6 0.8. 3.6. The outlet The velocity magnitude plots of the outlet plane have not been validated against experimental results.

960

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964 2

1.5

K 1 ⎯ UB Standard k - ω SST k - ω Standard k - ε RNG k - ε RSM

0.5

0

0

0.5

y ⎯ ⎯√Α

1

Fig. 31. Outlet: Mean velocity magnitude profiles normalised by bulk velocity, 10 l/min.

Nonetheless, they indicate that uniform or top-hat-shaped inlet conditions typically specified in lung airflow studies are incorrect. The flow patterns generated by the laryngeal jet are carried into the lower airways; lung studies should incorporate realistic inlet conditions in order to obtain physically realistic results. The velocity profiles of Fig. 31, as well as the three-dimensional velocity vectors plot of Fig. 32 show that, even for the lowest physiologically accurate flow rate of 10 l/min, unique inlet conditions (rather than uniform, parabolic, or turbulent velocity profiles) should be implemented for lower airways studies.

3.7. The three-dimensional flow Although the central sagittal plane velocity profiles are excellent arbiters of the accuracy of numerical results, they provide insufficient information to fully characterise the flow. Evidence provided by the smoke visualisations of [1] that the flow is highly three-dimensional in nature is supported the RANS modelled results. The Dean-type instabilities near station 3, as shown in Figs. 33, and 34 illustrate such phenomena. These corner-vortices are carried into and modified by the oro-pharynx, as shown in Fig. 35. This is immediately upstream of the backward facing step of the oro-pharynx, which is in turn upstream of station 5. Like the secondary vortices of the oral cavity in Fig. 33, flow structures in the pharynx, larynx, and trachea are effectively educed in iso-surface plots of negative second eigenvalue [36]. As labelled ‘a’ in Fig. 36, the recirculation bubble of the pharyngeal BFS creates flow structures that are carried the length of the pharynx to interact with the

Fig. 32. Outlet, coronal plane: Velocity vectors coloured by magnitude, RSM 10 l/min.

flow around the epiglottis. Also of interest are the recirculation bubble that is seen anterior to the epiglottis (b), and the shear layer and structures formed as the flow navigates the curved leading edge of the epiglottis (c). The three-dimensional nature of the flow at station 7, in the round-pipe trachea, is shown in Fig. 32. In this figure, the jet flow of the anterior portion of the trachea and the flow reversal of the posterior portion are seen. The impingement of the laryngeal jet on the anterior wall of

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

961

Fig. 33. Oro-pharynx: Instantaneous axial flow visualisation showing secondary flows, near station 3, taken from [24], 60 l/min.

Fig. 34. Oro-pharynx: iso-surface plot of k2 eigenvalue, coloured by vorticity magnitude, k–x, 10 l/min.

Fig. 35. Oro-pharynx, transverse plane. Top: instantaneous axial flow visualisation showing 2° flows, near station 4, 60 l/min from [1].

the trachea causes circumferential flow too, which results in repeating secondary vortices being formed, as shown in Fig. 37. These flow structures are carried, to a greater or lesser extent depending on flow rate, to the outlet of the trachea. This is critical information for those who study flow in the lower airways.

These iso-surface plots of negative second eigenvalue, from the k–x model results, provide significant insight into the nature of ETA flow. Analysis of complete time-dependent data sets is underway. Given the improved predictions of the mean velocity field, it is expected that the timedependent nature of the structures will provide significant

962

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

Fig. 36. Iso-surface of negative second eigenvalue, coloured by U-velocity, k–x model, 10 l/min.

insight into the effect of ETA morphology on the formation of secondary flows downstream. 4. Conclusion Given the tremendous variation in the quality of computational predictions across the measurement stations, a single defining or general statement about any of the evaluated turbulence models cannot be made. The relatively good performance of the k–x model at some, but not all, stations is circumspect; it is disconcerting that poor flow predictions at mid-stream locations had little effect on the outcome of downstream simulation results. Other turbulence models (the k–, the RNG k–x, and the Reynolds stress model) are not as consistently accurate in modelling of this internal complex flow. Nonetheless, the 10-fold increase in mesh refinement resulted in a greatly improved capacity to predict the flow phenomena in the ETA, as compared to earlier CFD simulations. To the extent that they are understood, the general flow features are well predicted by the k–x model. For wall bounded flows, it is normally expected that the k–x will outperform others because of its treatment of the viscous

Fig. 37. Iso-surface of negative second eigenvalue coloured by U-velocity; from k–x model, 10 l/min.

near-wall region and its handling of streamwise pressure gradients [22]. Beyond this, however, there is insufficient knowledge upon which to base any performance analysis for any turbulence model used to calculate flow in this complex geometry; there exists no documentation on CFD calculations for any of the combinations of canonical flow geometries (such as a curved leading edge, a BFS in a bent pipe, a variable height BFS, etc.) found in the ETA. Similarly, the physical limitations concomitant with experimental methods pose further hindrance to the ability to precisely describe the details of the flow. The highly three-dimensional, and probable time varying nature of the flow, as discovered by numerical methods, interacts with measurement probes used to measure the flow; additionally, the measurements that were taken in the central sagittal plane provide no information about the lateral secondary flows. It is vital that these flow phenomena be fully understood in order to confidently validate CFD models. Until the base flow can be modelled unquestioningly accurately, it would be imprudent to undertake efforts to model increasingly complex flow situations (such as cases with modified boundary conditions, or multi-phase flows). Should such studies be attempted, it would be impossible to identify the causes of inaccuracies in calculated solutions

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

– whether they be associated with the turbulence models, with the multi-phase models, with the wall treatment methods, and/or with other issues. Careful use of CFD methods, in conjunction with experimental methods, has permitted development of an improved understanding of the nature of ETA flow that would otherwise not be possible (with either method-type alone). The predicted main flow features – such as recirculation regions, secondary flow patterns, jets and impinging flows – are in accord with experiment. It is now possible to begin categorisation of the flow patterns associated with each portion of the ETA. At this point, however, only preliminary categorisation and understanding of these patterns is possible. There is no question as to the complexity of the flow inside the oral cavity (resulting from the unusual inlet and compound curved expansion), or of the flow around the epiglottis (a three-dimensional curved sharp leading edge), for example. Finally, it is clear that the complex flows developed in the ETA are carried into the trachea and to the lower airway. Knowing that tracheal flow is three-dimensional, those who are conducting studies of flow into the lungs should take care to apply realistic inlet conditions in their models. Clearly, the use of ‘standard’ inlet conditions (such as uniform, parabolic, or turbulent developed velocity profiles) would be inaccurate. Acknowledgements Support provided by the Natural Science and Engineering Research Council of Canada, and by the Queen’s University Graduate Awards is gratefully acknowledged. The authors thank the High Performance Computing Virtual Laboratory, Queen’s University site, and C3.ca for computational support. And importantly, the authors are grateful for the dedicated efforts of Emily Sheridan who produced some of the figures. References [1] Johnstone A, Uddin M, Pollard A, Heenan A, Finlay WH. The flow inside an idealised form of the human extra-thoracic airway. Exp Fluids 2004;37(5):673–89. [2] Heenan AF, Matida E, Pollard A, Finlay WH. Experimental measurements and computational modeling of the flow field in an idealized extra-thoracic airway. Exp Fluids 2003;35:70–84. [3] Gonda I. Targeting by deposition. In: Hickey AJ, editor. Pharmaceutical inhalation aerosol technology, vol. 54. New York: Marcel Dekker; 1992. p. 61–82. [4] Graham DI, James PW. Turbulent dispersion of particles using eddy interaction models. Int J Multiphase Flow 1996;22(1):157–75. [5] Hutchinson P, Hewitt GF, Dukler AE. Deposition of liquid or solid dispersions from turbulent gas streams: a stochastic model. Chem Eng Sci 1997;26(3):419–39. [6] Zhang Z, Kleinstreuer C. Low-Reynolds-number turbulent flows in locally constricted conduits: a comparison study. AIAA J 2003;41(5): 831–40. [7] Zhang Z, Kleinsteuer C. Species heat and mass transfer in a human upper airway model. Int J Heat Mass Transfer 2003;45:4755–68.

963

[8] Kleinstreuer C, Zhang Z. Laminar-to-turbulent fluid-particle flows in a human airway model. Int J Multiphase Flow 2003;29:271–89. [9] Zhang Z, Kleinstreuer C, Donohue J, Kim C. Comparison of microand nano-size particle depositions in a human upper airway model. Aerosol Sci 2005;36:211–33. [10] Matida E, Finlay W, Lange C, Grgic B. Improved numerical simulation of aerosol deposition in an idealised mouth–throat. Aerosol Sci 2004;35:1–19. [11] DeHaan WH, Finlay WH. Predicting extrathoracic deposition from dry powder inhalers. J Aerosol Sci 2004;35(3):309–31. [12] Heyerichs K, Pollard A. Heat transfer in separated and impinging turbulent flows. Int J Heat Mass Transfer 1996;39(12):2385–400. [13] Pollard A, Martinuzzi R. Comparative study of turbulence models in predicting turbulent pipe flow. Part II: Algebraic stress and k– models. AIAA J 1989;27(1):29–36. [14] Pollard A, Martinuzzi R. Comparative study of turbulence models in predicting turbulent pipe flow. Part II: Reynolds stress and k– models. AIAA J 1989;27(12):1714–21. [15] Poroseva S, Iaccarino G. Simulating separated flows using the k– model. In: Annual research briefs. Centre for Turbulence Research; 2001. p. 375–83. [16] Heenan A, Grgic B, Pollard A, Finlay W, Burnell PKP. An investigation of the relationship between the flow field and regional deposition in realistic extra-thoracic airways. J Aerosol Sci 2004;35(8):1013–23. [17] Stahlhofen W, Rudolf G, James AC. Intercomparison of experimental regional aerosol deposition data. J Aerosol Med 1989;2:285–308. [18] Stapleton KW, Guentsch E, Hoskinson MK, Finlay WH. the suitability of k– turbulence modeling for aerosol deposition in the mouth and throat: a comparison with experiment. J Aerosol Sci 2000;31(6):739–49. [19] Rhie CM, Chow WL. Numerical study of the turbulent flow past an airflow with trailing edge separation. AIAA J 1983;21:1525–32. [20] Launder BE, Spalding DB. Mathematical models of turbulence. New York: Academic Press; 1972. [21] Launder BE, Spalding DB. The numerical computation of turbulent flows. Comp Methods App Mech Eng 1974;3:269–89. [22] Wilcox DC. Turbulence modelling for CFD. 2nd ed. DCW Industries; 1993. [23] Chang HK, Masry OAE. A model study of flow dynamics in human central airways. Part 1: Axial velocity profiles. Resp Physiol 1982;49: 75–95. [24] Johnstone A. Hot wire measurements in an oropharyngeal pathway. MSc thesis. Kingston (ON), Canada: Queen’s University; November 2000. [25] Womersley JR. Method for the calculation of velocity rate of flow and viscous drag in arteries when the pressure gradient is known. J Physiol Lond 1955;127:553–63. [26] Finlay W. The mechanics of inhaled pharmaceutical aerosols. New York: Academic Press; 2001. [27] King M, Chang HK, Weber ME. Resistance of mucus-lined tubes to steady and oscillatory airflow. J Appl Physiol 1982;52(5):1172–6. [28] Hughes JMB, Hoppin FG, Mead J. Effect of lung inflation on bronchial length and diameter in excised lungs. J Appl Physiol 1972;32(1):25–35. [29] Jang YJ, Leschziner MA, Abe K, Temmerman L. Investigation of anisotropy-resolving turbulence models by reference to highlyresolved LES data for separated flow. Flow, Turb, Comb 2002;69:161–203. [30] Leschziner MA. Turbulence modelling for separated flows with anisotropy-resolving closures. Phil Trans R Soc Lond A 2000;358: 3247–77. [31] George WK, Abrahamsson H, Eriksson J, Karlsson RI, Lo¨fdahl L, Wosnik M. A similarity theory for the turbulent plane wall jet without external stream. J Fluid Mech 2000;425:367–411. [32] Lighthill MJ. Physiological fluid dynamics: a survey. J Fluid Mech 1972;52:475–97.

964

C.G. Ball et al. / Computers & Fluids 37 (2008) 943–964

[33] Corcoran TE, Chigier N. Characterization of the laryngeal jet using phase Doppler interferometry. J Aerosol Med 2000;13(2): 125–37. [34] Renotte C, Bouffious V, Wilquem F. Numerical 3d analysis of oscillatory flow in the time-varying laryngeal channel. J Biomech 2000;33:1637–77.

[35] Pedley TJ, Schroter RC, Sudlow MF. Pulmonary fluid dynamics. Ann Rev Fluid Mech 1977;9:74–229. [36] Jeong J, Hussain F. On the identification of a vortex. J Fluid Mech 1995;285:69–94. [37] Grant JCB. Grant’s atlas of anatomy. 5th ed. Baltimore: The Williams and Wilkins Co.; 1962.