Annals of BiomedicalEngineering,Vol. 20, pp. 573-581, 1992 Printedin the USA. All rights reserved.
0090-6964/92 $5.00 + .00 Copyright 9 1992Pergamon Press Ltd.
Numerical Simulation of Aerosol Particle Transport by Oscillating Flow in Respiratory Airways J a m e s K. Briant,* Duane D. F r a n k , * A n t h o n y and L. Loren Eylert
C. James,*
*Health Physics Department 1-Energy Sciences Department Pacific Northwest Laboratory Battelle Boulevard, Richland, WA
(Received3/7/91) Particle transport by oscillat&g flow & a tapered channel or in a tapered tube was computed from the complete equations of motion. These geometries represent a simplified model of the divergent flow field o f the mammalian bronchial tree. The computed deformation profile of a line of particles, transported by the oscillatory motion, was compared with prior experimental results and analytical calculations. All three methods agree that there is transport in the divergent direction of the tube by an axial stream of steady drift #n the core for moderately high frequency of oscillation (Womersley parameter in the range of I to 10). Bidirectional flow is established by an annular stream in the convergent direction, with no net flow on integral cycles of the oscillating fluid. A t higher frequency, however, the steady stream transforms to a different shape in the tapered tube, with transport in the divergent direction nearer the walls of the tube, rather than in the core. Transport by the continuing streams with oscillatory ventilation of the respiratory tract should deliver medicinal aerosols of low intrinsic particle mobility to the peripheral regions of the lungs. Keywords-High-frequency oscillatory ventilation, Mammalian bronchial tree, Computational fluid mechanics. INTRODUCTION T r a n s p o r t o f particles by oscillating m o t i o n in a t h r e e - d i m e n s i o n a l t a p e r e d t u b e is a s i m p l e r e p r e s e n t a t i o n o f p a r t i c l e t r a n s p o r t in the m a m m a l i a n lung. T h e f l o w field in a i r w a y s o f the b r o n c h i a l tree diverges as t h e a i r w a y s b i f u r c a t e in the m o r e p e r i p h eral regions o f the lungs, owing to the e x p a n s i o n o f t o t a l a i r w a y cross-sectional area f r o m 3 c m 2 in the t r a c h e a to m o r e t h a n 10,000 c m 2 in t h e a l v e o l a r d u c t s (8). By t h e 10th g e n e r a t i o n o f bifurcation, the cross-sectional a r e a o f a " p a r e n t " a i r w a y e x p a n d s
Acknowledgment-This work was supported by the U.S. Department of Energy under Contract DEAC06-76RLO 1830. The work of D. D. Frank was supported by the Northwest College and University Association for Science (Washington State University) under Grant DE-FG06-89ER-75522 with the U.S. Department of Energy. Address correspondence to James K. Briant, Health Physics Department, Pacific Northwest Laboratory, Battelle Boulevard, Richland, WA 99352. 573
574
J.K. Briant et al.
approximately 50% (from the middle of the parent airway to the sum of the two filial airways) and by the 20th generation, 800/0 (6). The expansion of cross-sectional area at a bifurcation is very gradual and smooth, starting about midway along the parent airway (7). Two basic flow phenomena occur in an airway bifurcation: gradual divergence of the flow field, and gradual change of direction. Modeling an airway bifurcation by using cylindrical tubes, which Simply intersect at the bifurcation, creates both flowfield divergence and change of direction, but in a rather abrupt manner. A standard airway bifurcation geometry with gradual divergence has not been established. Moreover, a mathematical model of the divergent flow-field for a viscous fluid in a bifurcation would be formidable. As an alternative, a simple tapered tube as a model of the bronchial airway greatly simplifies the mathematics of flow description by involving axial symmetry, which eliminates the directional component of the flow field in a bifurcation. The tapered tube model can be used to identify the role of flow-field divergence as the airway network of the bronchial tree expands in the periphery of the lungs. Analysis of the fluid mechanics in a tapered channel has been verified by experimental observation, and now by computer simulation. Using a regular perturbation method, Grotberg (4) derived analytical equations for oscillatory motion in a twodimensional tapered channel. A lubrication theory was involved in the derived description of oscillatory motion of a fixed stroke volume, with both fluid viscosity and inertia included. Pressure was assumed to be constant across the channel with a small angle of divergence. The slope of the tapered wall (e) was determined by the ratio of the total change in mean radius to total length of the corresponding airway. Grotberg chose the value e = 0.05 (2.86 ~ as being typical of airways. Continuing bidirectional axial streams (steady drifts) are predicted during oscillation of the fluid volume in the channel, with no net flow on integral cycles of oscillation. The continuing streams are in the divergent direction in the center of the channel (core of a tube), and in the convergent direction near the walls. The steady bidirectional drift predicted by Grotberg's theory (4) depends on the Womersley parameter (o0, up to a value where a double boundary layer develops, when the steady drift velocity profile changes shape. The Womersley parameter is the dimensionless ratio b(27rf/~,) ~/2, where b is the channel half-width (or radius of a tube), f is the oscillation frequency, and ~ is the kinematic viscosity. The bidirectional stream in a three-dimensional tapered tube can be expected to qualitatively resemble that in a two-dimensional tapered channel, but with a different dependence on the Womersley parameter. Gaver and Grotberg (2) experimentally produced pictures of transport in a tapered channel that was effectively two-dimensional, corresponding to Grotberg's analysis (4). Mixtures of glycerine and water were oscillated in the model channel for several values of the Womersley parameter. It was found that the deformation profile of a line of dye across the channel (perpendicular to the plane of symmetry) had been accurately described by the theory of Grotberg (4). Gaver and Grotberg also observed a start-up transient in the deformation profile, which eventually gave way to the steady bidirectional drift with continued oscillation. Grotberg's theory describes the bidirectional stream that develops beyond the influence of the start-up transient. At higher Womersley parameter (or > 10), development of the double boundary layer produced a major effect on the shape of the deformation profile, also as predicted by Grotberg.
Aerosol Transport by Oscillating Flow
575
Godleski and Grotberg (3) made a theoretical analysis of transport by oscillatory motion in a three-dimensional tapered tube. Their analysis also included the diffusion parameter to model transport of gases. The steady drift velocity at the center-line was calculated to be about 5 times greater in the tapered tube than in a corresponding tapered channel for o~ in the range of 0.5 to 5. Motion was focused along the centerline of the three-dimensional tube and amplified by the geometry, a more realistic model of the respiratory tract airways than the two-dimensional channel. In particular, transport in the tapered tube represents transport in the lung airways for the short-stroke oscillatory motion of high-frequency ventilation, where the local stroke length is small relative to the dimensions of a bifurcation. Transport of 0.5-/~m aerosol particles in a tapered tube is simulated by computation of nondiffusional motion. This particle size has minimal intrinsic mobility, owing to the transition from thermal mobility to inertial mobility at around 0.5-/~m particle diameter, and can be considered a tracer of convection. Such transport was observed by Briant and Lippmann (I) with oscillation of various gases in a hollow cast of the tracheobronchial tree. Aerosol particles were rapidly transported from the trachea to small terminal airways by high-frequency (10- and 21-Hz) oscillation of volumes with one-fifth that of the cast. Gases of higher kinematic viscosity gave faster, more efficient transport of the particles, demonstrating transport by axial streaming during high-frequency oscillatory ventilation. The present investigation validates computer simulations of transport profiles which match the two-dimensional transport experiments of Gaver and Grotberg (2), and is extended to visualize transport in a three-dimensional tapered tube. This study compares the profile of deformation of a straight line of simulated particles across the two-dimensional tapered channel to that of deformation of a flat plane across a three-dimensional tapered tube for corresponding dimensions. The computer simulations include start-up transients which gradually become insignificant with continuous oscillation. After a number of oscillations, the shapes of deformation profiles agree with the analytical calculations derived by the methods of Grotberg (4) and Godleski and Grotberg (3). Visualizing the fluid mechanics of oscillatory motion in a tapered tube by analytical calculation or computational simulation gives us better understanding of transport through airways of the lung. METHODS
The computer code TEMPEST 1 was used to perform the numerical simulation of flow fields in tapered channels and tubes. TEMPEST uses a time-marching finite-difference procedure to solve equations of motion based on the Navier-Stokes equations for the flow field of interest. The governing equations of motion, including continuity, momentum, and transport are solved sequentially. The governing equations in generalized curvilinear coordinates are as follows:
Equation o f Continuity 0 (4~v;) OTf
--~
= o 9
1Transient Energy, Momentum,and Pressure Equation Solution in Three dimensions.Proprietary software of Battelle,PacificNorthwestLaboratories.
576
J.K. Briant et al. Equation o f M o m e n t u m (Covariant Form)
ave
-
-
1 -~-
m
Ot = --
0
a~lj (~/gViVJ) -- [ i j , m ] V m V i.j
m
1 aP
1
m + -no 0~1' ~
O
(
- - (~,/~ve/) n o~J ~ j
i )
vegn +
P gi --
Po
. -- Di
,
Equation o f Transport (Scalar) adp 1 a a-7 + --~ ---~ a,7 ( ~ v , )
-
x/gl O~l kO ( 4 ~ g j % ~ O ~ ) + S , .
where, for indices i, j, k, m, n: 17i ~--- coordinate transform space g = Jacobian V i = contravariant velocity component Vi = covariant velocity component t = time [ij, m] = Christoffel symbol of the first kind p = density (with po being at time 0) P = pressure p = kinematic viscosity e / = strain tensor In
ij
/ = Christoffel symbol of the second kind ,i
Errg
g~ = Dg = q~ = gU = % = S, =
strain tensor gravitational vector flow drag term scalar transported metric tensor scalar diffusion coefficient scalar source.
In addition, options are available in T E M P E S T to solve turbulence model equations and treat non-Newtonian fluids. The solution procedure used in T E M P E S T is a semi-implicit time-marching scheme. At each time step, m o m e n t u m equations are solved explicitly, and pressure implicitly. Scalar transport equations are solved explicitly or implicitly, depending upon stability criteria. Computation grid structure is required in the form of a boundary-fitted curvilinear coordinate system. Fluid properties, boundary conditions, and simulation control are user-supplied. For the present computations, temperature and pressure were fixed at 20~ and 1 atm. The tapered channel used for the present computation had the same configuration and dimensions as the one used experimentally by Gaver and Grotberg (2). The tapered section o f the channel was 20 cm long and 60 cm high, expanding f r o m 3 cm width at one end to 5 cm at the other end, with straight sections extending from each end. The slope of the walls of the tapered section was e = 0.05 (2.86~ A line of par-
Aerosol Transport by Oscillating Flow
577
ticles was positioned across the middle of the channel where the channel half-width (b) is 2 cm. These particles were given no inertia or diffusibility, and hence served simply as convective tracers. Positions of the particles were graphically observed at integral cycles of oscillation, with the piston at midstroke. Piston motion was sinusoidal in all of the models. Although particle positions could be tracked in all computation cells, only those that originated in the layer of cells across the middle of the taper are discussed in this article. The tapered tube had the same dimensions as the tapered channel, with the radius of the tube section equal to the half-width (b) of the channel section. Thus, a diametrical section of the tapered tube looked just like a section of the tapered channel. Piston stroke amplitude was the same for the equivalent channel and tube comparisons. The dimensionless amplitude was ,4 = d/b, where d is the characteristic stroke length (stroke volume divided by cross-sectional area) at the reference position at half-width b = 2 cm. RESULTS The computed particle transport in the tapered channel closely matched the shape of that experimentally observed by Gaver and Grotberg (2), as shown in Fig. 1. Total number of oscillations could be adjusted to produce the same length of deformation shown in the photograph by Gaver and Grotberg (2, see their figure 2a). The TEMPEST simulation accommodated computation of the start-up transient that occurred in the experimental model. Figure 2 compares the TEMPEST computation of the start-up transient with the measurement by Gaver and Grotberg (2). In Fig. 2, the dimensionless amplitude was ,4 = 2.38, while all other calculations and computations reported here used ,4 = 1.59. W 1" is the total length of the deformation profile in the axial direction, and increases with continued oscillation (the asterisk denotes total cumulation). Computation of W 1" for the first cycle is strongly influenced by the
.........................................
TEMP~
I
/ Gaverand Grotberg FIGURE 1. Computed deformation profile includes the start-up transient for start at mid-stroke in the divergent direction, with A = 1 . 5 9 , e = 0 . 0 5 , and a = 3 . 8 f o r a total of six complete strokes in the two-dimensional channel (similar t o the experimental deformation profile shown by Gaver and Grotherg [2]; see their figure 2a).
578
J.K. Briant et al.
7
9 .0 9
9 Gaver & Grotberg data o-
6
o TEMPEST computation
.0" .0 ~ **
5 ~
E ._
D"
.d
4
W l * = -1.5 + 0.028t
0. 9
3
o.0" 9~ 9 9149
2
.o . 9 1 4~176 9
Q."
1 0 -1
i 9 .-'"
",.."i o ."
t
i
i
i
i
100
i
i 200
i
i
=
i
t 300
Time (s)
99 O
FIGURE 2. Start-up of transport by sinusoidal oscillations, as observed by Gaver and Grotberg (2) and computed by TEMPEST, is plotted for distance along the center line of a two-dimensional channel, w i t h A = 2 . 3 8 , ( = 0 . 0 5 , a n d a = 14.
initial acceleration of the piston at the time of start-up, which can only be estimated for computation of the start-up transient. Particle transport in a tapered tube, and in a tapered channel of corresponding dimensions, is shown in Fig. 3 for deformation of a straight line by one cycle of oscillation during continuous oscillation for both analytical calculation and numerical computation. Axial transport by the steady drift for each oscillation is given as W1 S, defined as the length of deformation for each oscillation during continuous oscillation. This is plotted for the steady drift as a fraction of tube radius or channel halfwidth9 Calculated transport in the two-dimensional channel is nearly indistinguishable from the computed transport, which takes the start-up transient into account. Calculated transport in the three-dimensional channel is a shape similar to the computed transport, but of less magnitude. The difference in magnitude might be caused by the start-up transient, which has not yet been investigated for the three-dimensional case. At higher values of Womersley parameter (o~ > 10), a double boundary layer develops and changes the form of the steady drift. Gaver and Grotberg (2) observed this in the two-dimensional tapered channel and found it to be consistent with the prediction from the analytical calculations given by Grotberg (4). The double boundary layer is a complex challenge that TEMPEST handled well. Figure 4 demonstrates transport in the same channel shown in Fig. 1, but at a higher value of Womersley parameter, or = 17.4; Figs. 1 and 4 both show the total transport for six oscillations to allow direct comparison for two different values of a. The difference between Fig. 4 and the photograph given by Gaver and Grotberg (2, see their figure 2b) could be caused by differences in the start-up conditions, which were not defined in their
Aerosol Transport by Oscillating Flow 0.5
579
m
~~,~,.
0.4 0.3
~
E
'
M
P
E
3D', Godle'ski& G'rotberg' S
'
'
T
0.2 E
(,3 e'}
~
0.1 0.0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
'
.
.
.
2D, Grotberg
.
-0.1 -0.2 0.0
I
I
0.2
I
I
0.4
I
r (cm)
I
0.6
I
I
0.8
I
1.0
Comparison of analytical solution and numerical computation by TEMPEST is made for and three-dimensional cases of similar dimensions for each oscillation during steady drift (WlS). Fractional distance from the center-line axis is r for either the channel or the tube (which coincides with distance in centimeters for these models). Analytical solutions are for a = 5 . 0 ; numerical computations, for a = 3 , 8 , The three-dimensional case does not account for start-up transient. FIGURE 3.
two-
article. As in Fig. 1, the length of the deformation profile could be adjusted by the total number of oscillations. Most importantly, TEMPEST does exhibit the double boundary layer in the same range of Womersley parameter predicted by Grotberg (4) and observed by Gaver and Grotberg (2). Computations by TEMPEST at higher Womersley parameter values in the three-dimensional tapered tube, however, do not properly converge in some regions of the model. We hope eventually to resolve this problem through further investigation. DISCUSSION
The start-up transient is remarkably varied for different positions in the oscillation cycle and direction of start of the piston. When started in opposite directions from the same position, the deformation profiles are substantially different at the end of the first cycle. This demonstrates the importance of inertia along with a lubrication theory involving an unsteady flow field. If the fluid inertia is assumed to be negligible relative to fluid viscosity, as in Stokes flow, then the deformation profiles are expected to be the same, regardless of start-up conditions. Therefore, an assumption of laminar parallel flow is not appropriate for the tapered airway where flow is typically unsteady, as in the tapered tube model, Inertia of the fluid is not negligible, except at extremely low velocities in steady flow (5). When flow is unsteady, inertial effects result in quasi-steady streaming during oscillatory motion, even with Womersley parameter considerably less than 1 (3).
580
J.K. Briant et al.
~
'
TEMPEST
Gaver and Grotberg
|
|
FIGURE 4. Computed deformation profile includes the start-up transient for start at mid-stroke in the divergent direction, with A = 1 . 5 9 , ~ = 0 . 0 5 , and a = 1 7 . 4 f o r a t o t a l o f six complete strokes in the two-dimensional channel (similar to Gavar and Grotberg [2]; see their figure 2b). Although the start-up transients may have been different, both exhibit the double boundary layer.
Transport in the three-dimensional tube shows a substantially enhanced axial core stream, compared to the two-dimensional channel, as shown in Fig. 3. Maximum amplitude of transport has not been calculated as a function of Womersley parameter. Godleski and Grotberg (3) calculated that gas transport increases with Womersley parameter to value 12 and beyond. At higher Womersley parameter values (o~ > 10), transport in the divergent direction is in the annular region of the airway, while transport in the core is in the convergent direction. The transition to the double boundary layer form of the transport (deformation) profile occurs at a Womersley parameter value of about 10 in the two-dimensional channel. The difference between the two forms can be seen by comparing Figs. 1 and 4, which differ only in viscosity of the fluid. The steady axial streaming during high-frequency ventilation should be present in the mammalian bronchial tree, which has a flow-field divergence similar to the tapered tube. Molecules of an insoluble gas should be transported from the most central airway (trachea) to the most peripheral (alveolar ducts) by the steady bidirectional stream through a wide range of Womersley parameter values. At a high Womersley parameter value (a > 10), gas molecules would be transported by the annular stream in the divergent direction near the airway walls. Aerosol particles in this region would be more likely to interact with the airway surface by deposition on the walls. Thus, transport of aerosol particles should be more efficient at a lower Womersley parameter value, when the particles would be carried in the divergent direction in the axial core stream isolated from the airway walls. Even when the Womersley parameter is as high as 30 in the trachea, it is less than 10 through much of the bronchial tree. It appears that high-frequency oscillatory ventilation could be used to efficiently deliver medicinal aerosol particles to the peripheral regions of the lungs. Bifurcation geometry undoubtedly influences the structure of the axial bidirectional stream in the bifurcating network of the bronchial tree. Short-stroke high-frequency oscillation produces motion smaller than the characteristic dimension of a bifurca-
Aerosol Transport by Oscillating Flow
581
t i o n , which leaves a i r w a y t a p e r as the d o m i n a n t f a c t o r ; f o r such oscillations, t h e n , d i r e c t i o n a l g e o m e t r y has m i n i m a l influence o n t r a n s p o r t . T h e l o n g e r - s t r o k e , l o w e r frequency oscillation o f n o r m a l ventilation involves lengths c o n s i d e r a b l y greater t h a n a n y d i m e n s i o n o f a b i f u r c a t i o n , thus a d d i n g m o r e d i r e c t i o n a l d e p e n d e n c e ; the t r a n s p o r t p r o f i l e is skewed f r o m the a i r w a y center line in d i f f e r e n t ways with e a c h direct i o n o f flow t h r o u g h a bifurcation. Thus, the p e n e t r a t i o n o f gas a n d a e r o s o l particles to t h e p e r i p h e r a l airways o f the lung b y n o r m a l v e n t i l a t i o n does i n v o l v e b i f u r c a t i n g f l o w - f i e l d g e o m e t r y ; to c o m p l e t e l y u n d e r s t a n d this p r o c e s s , the simple t a p e r e d t u b e m o d e l is p r o b a b l y n o t a d e q u a t e . E v e n t u a l l y , T E M P E S T s h o u l d be a b l e to c o m p u t e t r a n s p o r t in b i f u r c a t i n g flow fields t h a t b e t t e r m o d e l a i r w a y s o f t h e lungs.
REFERENCES 1. Briant, J.K.; Lippmann, M. Particle transport through a hollow canine airway cast by high-frequency oscillatory ventilation. Exp. Lung Res. 18:385-407; 1992. 2. Gaver III, D.P.; Grotberg, J.B. An experimental investigation of oscillating flow in a tapered channel. J. Fluid Mech. 172:47-61; 1986. 3. Godleski, D.A.; Grotberg, J.B. Convection-diffusion interaction for oscillatory flow in a tapered tube. J. Biomech. Eng. 110:283-291; 1988. 4. Grotberg, J.B. Volume-cycled oscillatory flow in a tapered channel. J. Fluid Mech. 141:249-264; 1984. 5. Happel, J.; Brenner, H. Low Reynolds number hydrodynamics, with special applications to particulate media. The Hague, The Netherlands: Martinus Nijhoff Publishers; 1983. 6. Horsfield, K.; Cumming, G. Angles of branching and diameters of branches in the human bronchial tree. Bull. Math. Biophys. 29:245-259; 1967. 7. Schreck, R.M. Laminar flow through bifurcations with applications to the human lung. Ph.D. Dissertation. Evanston, IL: Northwestern University; 1972. 8. Weibel, E.R. Morphometry of the human lung. New York: Academic Press; 1963.
0090-6964/92 $5.00 + .00 Copyright 9 1992Pergamon Press Ltd.
Numerical Simulation of Aerosol Particle Transport by Oscillating Flow in Respiratory Airways J a m e s K. Briant,* Duane D. F r a n k , * A n t h o n y and L. Loren Eylert
C. James,*
*Health Physics Department 1-Energy Sciences Department Pacific Northwest Laboratory Battelle Boulevard, Richland, WA
(Received3/7/91) Particle transport by oscillat&g flow & a tapered channel or in a tapered tube was computed from the complete equations of motion. These geometries represent a simplified model of the divergent flow field o f the mammalian bronchial tree. The computed deformation profile of a line of particles, transported by the oscillatory motion, was compared with prior experimental results and analytical calculations. All three methods agree that there is transport in the divergent direction of the tube by an axial stream of steady drift #n the core for moderately high frequency of oscillation (Womersley parameter in the range of I to 10). Bidirectional flow is established by an annular stream in the convergent direction, with no net flow on integral cycles of the oscillating fluid. A t higher frequency, however, the steady stream transforms to a different shape in the tapered tube, with transport in the divergent direction nearer the walls of the tube, rather than in the core. Transport by the continuing streams with oscillatory ventilation of the respiratory tract should deliver medicinal aerosols of low intrinsic particle mobility to the peripheral regions of the lungs. Keywords-High-frequency oscillatory ventilation, Mammalian bronchial tree, Computational fluid mechanics. INTRODUCTION T r a n s p o r t o f particles by oscillating m o t i o n in a t h r e e - d i m e n s i o n a l t a p e r e d t u b e is a s i m p l e r e p r e s e n t a t i o n o f p a r t i c l e t r a n s p o r t in the m a m m a l i a n lung. T h e f l o w field in a i r w a y s o f the b r o n c h i a l tree diverges as t h e a i r w a y s b i f u r c a t e in the m o r e p e r i p h eral regions o f the lungs, owing to the e x p a n s i o n o f t o t a l a i r w a y cross-sectional area f r o m 3 c m 2 in the t r a c h e a to m o r e t h a n 10,000 c m 2 in t h e a l v e o l a r d u c t s (8). By t h e 10th g e n e r a t i o n o f bifurcation, the cross-sectional a r e a o f a " p a r e n t " a i r w a y e x p a n d s
Acknowledgment-This work was supported by the U.S. Department of Energy under Contract DEAC06-76RLO 1830. The work of D. D. Frank was supported by the Northwest College and University Association for Science (Washington State University) under Grant DE-FG06-89ER-75522 with the U.S. Department of Energy. Address correspondence to James K. Briant, Health Physics Department, Pacific Northwest Laboratory, Battelle Boulevard, Richland, WA 99352. 573
574
J.K. Briant et al.
approximately 50% (from the middle of the parent airway to the sum of the two filial airways) and by the 20th generation, 800/0 (6). The expansion of cross-sectional area at a bifurcation is very gradual and smooth, starting about midway along the parent airway (7). Two basic flow phenomena occur in an airway bifurcation: gradual divergence of the flow field, and gradual change of direction. Modeling an airway bifurcation by using cylindrical tubes, which Simply intersect at the bifurcation, creates both flowfield divergence and change of direction, but in a rather abrupt manner. A standard airway bifurcation geometry with gradual divergence has not been established. Moreover, a mathematical model of the divergent flow-field for a viscous fluid in a bifurcation would be formidable. As an alternative, a simple tapered tube as a model of the bronchial airway greatly simplifies the mathematics of flow description by involving axial symmetry, which eliminates the directional component of the flow field in a bifurcation. The tapered tube model can be used to identify the role of flow-field divergence as the airway network of the bronchial tree expands in the periphery of the lungs. Analysis of the fluid mechanics in a tapered channel has been verified by experimental observation, and now by computer simulation. Using a regular perturbation method, Grotberg (4) derived analytical equations for oscillatory motion in a twodimensional tapered channel. A lubrication theory was involved in the derived description of oscillatory motion of a fixed stroke volume, with both fluid viscosity and inertia included. Pressure was assumed to be constant across the channel with a small angle of divergence. The slope of the tapered wall (e) was determined by the ratio of the total change in mean radius to total length of the corresponding airway. Grotberg chose the value e = 0.05 (2.86 ~ as being typical of airways. Continuing bidirectional axial streams (steady drifts) are predicted during oscillation of the fluid volume in the channel, with no net flow on integral cycles of oscillation. The continuing streams are in the divergent direction in the center of the channel (core of a tube), and in the convergent direction near the walls. The steady bidirectional drift predicted by Grotberg's theory (4) depends on the Womersley parameter (o0, up to a value where a double boundary layer develops, when the steady drift velocity profile changes shape. The Womersley parameter is the dimensionless ratio b(27rf/~,) ~/2, where b is the channel half-width (or radius of a tube), f is the oscillation frequency, and ~ is the kinematic viscosity. The bidirectional stream in a three-dimensional tapered tube can be expected to qualitatively resemble that in a two-dimensional tapered channel, but with a different dependence on the Womersley parameter. Gaver and Grotberg (2) experimentally produced pictures of transport in a tapered channel that was effectively two-dimensional, corresponding to Grotberg's analysis (4). Mixtures of glycerine and water were oscillated in the model channel for several values of the Womersley parameter. It was found that the deformation profile of a line of dye across the channel (perpendicular to the plane of symmetry) had been accurately described by the theory of Grotberg (4). Gaver and Grotberg also observed a start-up transient in the deformation profile, which eventually gave way to the steady bidirectional drift with continued oscillation. Grotberg's theory describes the bidirectional stream that develops beyond the influence of the start-up transient. At higher Womersley parameter (or > 10), development of the double boundary layer produced a major effect on the shape of the deformation profile, also as predicted by Grotberg.
Aerosol Transport by Oscillating Flow
575
Godleski and Grotberg (3) made a theoretical analysis of transport by oscillatory motion in a three-dimensional tapered tube. Their analysis also included the diffusion parameter to model transport of gases. The steady drift velocity at the center-line was calculated to be about 5 times greater in the tapered tube than in a corresponding tapered channel for o~ in the range of 0.5 to 5. Motion was focused along the centerline of the three-dimensional tube and amplified by the geometry, a more realistic model of the respiratory tract airways than the two-dimensional channel. In particular, transport in the tapered tube represents transport in the lung airways for the short-stroke oscillatory motion of high-frequency ventilation, where the local stroke length is small relative to the dimensions of a bifurcation. Transport of 0.5-/~m aerosol particles in a tapered tube is simulated by computation of nondiffusional motion. This particle size has minimal intrinsic mobility, owing to the transition from thermal mobility to inertial mobility at around 0.5-/~m particle diameter, and can be considered a tracer of convection. Such transport was observed by Briant and Lippmann (I) with oscillation of various gases in a hollow cast of the tracheobronchial tree. Aerosol particles were rapidly transported from the trachea to small terminal airways by high-frequency (10- and 21-Hz) oscillation of volumes with one-fifth that of the cast. Gases of higher kinematic viscosity gave faster, more efficient transport of the particles, demonstrating transport by axial streaming during high-frequency oscillatory ventilation. The present investigation validates computer simulations of transport profiles which match the two-dimensional transport experiments of Gaver and Grotberg (2), and is extended to visualize transport in a three-dimensional tapered tube. This study compares the profile of deformation of a straight line of simulated particles across the two-dimensional tapered channel to that of deformation of a flat plane across a three-dimensional tapered tube for corresponding dimensions. The computer simulations include start-up transients which gradually become insignificant with continuous oscillation. After a number of oscillations, the shapes of deformation profiles agree with the analytical calculations derived by the methods of Grotberg (4) and Godleski and Grotberg (3). Visualizing the fluid mechanics of oscillatory motion in a tapered tube by analytical calculation or computational simulation gives us better understanding of transport through airways of the lung. METHODS
The computer code TEMPEST 1 was used to perform the numerical simulation of flow fields in tapered channels and tubes. TEMPEST uses a time-marching finite-difference procedure to solve equations of motion based on the Navier-Stokes equations for the flow field of interest. The governing equations of motion, including continuity, momentum, and transport are solved sequentially. The governing equations in generalized curvilinear coordinates are as follows:
Equation o f Continuity 0 (4~v;) OTf
--~
= o 9
1Transient Energy, Momentum,and Pressure Equation Solution in Three dimensions.Proprietary software of Battelle,PacificNorthwestLaboratories.
576
J.K. Briant et al. Equation o f M o m e n t u m (Covariant Form)
ave
-
-
1 -~-
m
Ot = --
0
a~lj (~/gViVJ) -- [ i j , m ] V m V i.j
m
1 aP
1
m + -no 0~1' ~
O
(
- - (~,/~ve/) n o~J ~ j
i )
vegn +
P gi --
Po
. -- Di
,
Equation o f Transport (Scalar) adp 1 a a-7 + --~ ---~ a,7 ( ~ v , )
-
x/gl O~l kO ( 4 ~ g j % ~ O ~ ) + S , .
where, for indices i, j, k, m, n: 17i ~--- coordinate transform space g = Jacobian V i = contravariant velocity component Vi = covariant velocity component t = time [ij, m] = Christoffel symbol of the first kind p = density (with po being at time 0) P = pressure p = kinematic viscosity e / = strain tensor In
ij
/ = Christoffel symbol of the second kind ,i
Errg
g~ = Dg = q~ = gU = % = S, =
strain tensor gravitational vector flow drag term scalar transported metric tensor scalar diffusion coefficient scalar source.
In addition, options are available in T E M P E S T to solve turbulence model equations and treat non-Newtonian fluids. The solution procedure used in T E M P E S T is a semi-implicit time-marching scheme. At each time step, m o m e n t u m equations are solved explicitly, and pressure implicitly. Scalar transport equations are solved explicitly or implicitly, depending upon stability criteria. Computation grid structure is required in the form of a boundary-fitted curvilinear coordinate system. Fluid properties, boundary conditions, and simulation control are user-supplied. For the present computations, temperature and pressure were fixed at 20~ and 1 atm. The tapered channel used for the present computation had the same configuration and dimensions as the one used experimentally by Gaver and Grotberg (2). The tapered section o f the channel was 20 cm long and 60 cm high, expanding f r o m 3 cm width at one end to 5 cm at the other end, with straight sections extending from each end. The slope of the walls of the tapered section was e = 0.05 (2.86~ A line of par-
Aerosol Transport by Oscillating Flow
577
ticles was positioned across the middle of the channel where the channel half-width (b) is 2 cm. These particles were given no inertia or diffusibility, and hence served simply as convective tracers. Positions of the particles were graphically observed at integral cycles of oscillation, with the piston at midstroke. Piston motion was sinusoidal in all of the models. Although particle positions could be tracked in all computation cells, only those that originated in the layer of cells across the middle of the taper are discussed in this article. The tapered tube had the same dimensions as the tapered channel, with the radius of the tube section equal to the half-width (b) of the channel section. Thus, a diametrical section of the tapered tube looked just like a section of the tapered channel. Piston stroke amplitude was the same for the equivalent channel and tube comparisons. The dimensionless amplitude was ,4 = d/b, where d is the characteristic stroke length (stroke volume divided by cross-sectional area) at the reference position at half-width b = 2 cm. RESULTS The computed particle transport in the tapered channel closely matched the shape of that experimentally observed by Gaver and Grotberg (2), as shown in Fig. 1. Total number of oscillations could be adjusted to produce the same length of deformation shown in the photograph by Gaver and Grotberg (2, see their figure 2a). The TEMPEST simulation accommodated computation of the start-up transient that occurred in the experimental model. Figure 2 compares the TEMPEST computation of the start-up transient with the measurement by Gaver and Grotberg (2). In Fig. 2, the dimensionless amplitude was ,4 = 2.38, while all other calculations and computations reported here used ,4 = 1.59. W 1" is the total length of the deformation profile in the axial direction, and increases with continued oscillation (the asterisk denotes total cumulation). Computation of W 1" for the first cycle is strongly influenced by the
.........................................
TEMP~
I
/ Gaverand Grotberg FIGURE 1. Computed deformation profile includes the start-up transient for start at mid-stroke in the divergent direction, with A = 1 . 5 9 , e = 0 . 0 5 , and a = 3 . 8 f o r a total of six complete strokes in the two-dimensional channel (similar t o the experimental deformation profile shown by Gaver and Grotherg [2]; see their figure 2a).
578
J.K. Briant et al.
7
9 .0 9
9 Gaver & Grotberg data o-
6
o TEMPEST computation
.0" .0 ~ **
5 ~
E ._
D"
.d
4
W l * = -1.5 + 0.028t
0. 9
3
o.0" 9~ 9 9149
2
.o . 9 1 4~176 9
Q."
1 0 -1
i 9 .-'"
",.."i o ."
t
i
i
i
i
100
i
i 200
i
i
=
i
t 300
Time (s)
99 O
FIGURE 2. Start-up of transport by sinusoidal oscillations, as observed by Gaver and Grotberg (2) and computed by TEMPEST, is plotted for distance along the center line of a two-dimensional channel, w i t h A = 2 . 3 8 , ( = 0 . 0 5 , a n d a = 14.
initial acceleration of the piston at the time of start-up, which can only be estimated for computation of the start-up transient. Particle transport in a tapered tube, and in a tapered channel of corresponding dimensions, is shown in Fig. 3 for deformation of a straight line by one cycle of oscillation during continuous oscillation for both analytical calculation and numerical computation. Axial transport by the steady drift for each oscillation is given as W1 S, defined as the length of deformation for each oscillation during continuous oscillation. This is plotted for the steady drift as a fraction of tube radius or channel halfwidth9 Calculated transport in the two-dimensional channel is nearly indistinguishable from the computed transport, which takes the start-up transient into account. Calculated transport in the three-dimensional channel is a shape similar to the computed transport, but of less magnitude. The difference in magnitude might be caused by the start-up transient, which has not yet been investigated for the three-dimensional case. At higher values of Womersley parameter (o~ > 10), a double boundary layer develops and changes the form of the steady drift. Gaver and Grotberg (2) observed this in the two-dimensional tapered channel and found it to be consistent with the prediction from the analytical calculations given by Grotberg (4). The double boundary layer is a complex challenge that TEMPEST handled well. Figure 4 demonstrates transport in the same channel shown in Fig. 1, but at a higher value of Womersley parameter, or = 17.4; Figs. 1 and 4 both show the total transport for six oscillations to allow direct comparison for two different values of a. The difference between Fig. 4 and the photograph given by Gaver and Grotberg (2, see their figure 2b) could be caused by differences in the start-up conditions, which were not defined in their
Aerosol Transport by Oscillating Flow 0.5
579
m
~~,~,.
0.4 0.3
~
E
'
M
P
E
3D', Godle'ski& G'rotberg' S
'
'
T
0.2 E
(,3 e'}
~
0.1 0.0
.
.
.
.
.
.
.
.
.
.
.
.
.
.
'
.
.
.
2D, Grotberg
.
-0.1 -0.2 0.0
I
I
0.2
I
I
0.4
I
r (cm)
I
0.6
I
I
0.8
I
1.0
Comparison of analytical solution and numerical computation by TEMPEST is made for and three-dimensional cases of similar dimensions for each oscillation during steady drift (WlS). Fractional distance from the center-line axis is r for either the channel or the tube (which coincides with distance in centimeters for these models). Analytical solutions are for a = 5 . 0 ; numerical computations, for a = 3 , 8 , The three-dimensional case does not account for start-up transient. FIGURE 3.
two-
article. As in Fig. 1, the length of the deformation profile could be adjusted by the total number of oscillations. Most importantly, TEMPEST does exhibit the double boundary layer in the same range of Womersley parameter predicted by Grotberg (4) and observed by Gaver and Grotberg (2). Computations by TEMPEST at higher Womersley parameter values in the three-dimensional tapered tube, however, do not properly converge in some regions of the model. We hope eventually to resolve this problem through further investigation. DISCUSSION
The start-up transient is remarkably varied for different positions in the oscillation cycle and direction of start of the piston. When started in opposite directions from the same position, the deformation profiles are substantially different at the end of the first cycle. This demonstrates the importance of inertia along with a lubrication theory involving an unsteady flow field. If the fluid inertia is assumed to be negligible relative to fluid viscosity, as in Stokes flow, then the deformation profiles are expected to be the same, regardless of start-up conditions. Therefore, an assumption of laminar parallel flow is not appropriate for the tapered airway where flow is typically unsteady, as in the tapered tube model, Inertia of the fluid is not negligible, except at extremely low velocities in steady flow (5). When flow is unsteady, inertial effects result in quasi-steady streaming during oscillatory motion, even with Womersley parameter considerably less than 1 (3).
580
J.K. Briant et al.
~
'
TEMPEST
Gaver and Grotberg
|
|
FIGURE 4. Computed deformation profile includes the start-up transient for start at mid-stroke in the divergent direction, with A = 1 . 5 9 , ~ = 0 . 0 5 , and a = 1 7 . 4 f o r a t o t a l o f six complete strokes in the two-dimensional channel (similar to Gavar and Grotberg [2]; see their figure 2b). Although the start-up transients may have been different, both exhibit the double boundary layer.
Transport in the three-dimensional tube shows a substantially enhanced axial core stream, compared to the two-dimensional channel, as shown in Fig. 3. Maximum amplitude of transport has not been calculated as a function of Womersley parameter. Godleski and Grotberg (3) calculated that gas transport increases with Womersley parameter to value 12 and beyond. At higher Womersley parameter values (o~ > 10), transport in the divergent direction is in the annular region of the airway, while transport in the core is in the convergent direction. The transition to the double boundary layer form of the transport (deformation) profile occurs at a Womersley parameter value of about 10 in the two-dimensional channel. The difference between the two forms can be seen by comparing Figs. 1 and 4, which differ only in viscosity of the fluid. The steady axial streaming during high-frequency ventilation should be present in the mammalian bronchial tree, which has a flow-field divergence similar to the tapered tube. Molecules of an insoluble gas should be transported from the most central airway (trachea) to the most peripheral (alveolar ducts) by the steady bidirectional stream through a wide range of Womersley parameter values. At a high Womersley parameter value (a > 10), gas molecules would be transported by the annular stream in the divergent direction near the airway walls. Aerosol particles in this region would be more likely to interact with the airway surface by deposition on the walls. Thus, transport of aerosol particles should be more efficient at a lower Womersley parameter value, when the particles would be carried in the divergent direction in the axial core stream isolated from the airway walls. Even when the Womersley parameter is as high as 30 in the trachea, it is less than 10 through much of the bronchial tree. It appears that high-frequency oscillatory ventilation could be used to efficiently deliver medicinal aerosol particles to the peripheral regions of the lungs. Bifurcation geometry undoubtedly influences the structure of the axial bidirectional stream in the bifurcating network of the bronchial tree. Short-stroke high-frequency oscillation produces motion smaller than the characteristic dimension of a bifurca-
Aerosol Transport by Oscillating Flow
581
t i o n , which leaves a i r w a y t a p e r as the d o m i n a n t f a c t o r ; f o r such oscillations, t h e n , d i r e c t i o n a l g e o m e t r y has m i n i m a l influence o n t r a n s p o r t . T h e l o n g e r - s t r o k e , l o w e r frequency oscillation o f n o r m a l ventilation involves lengths c o n s i d e r a b l y greater t h a n a n y d i m e n s i o n o f a b i f u r c a t i o n , thus a d d i n g m o r e d i r e c t i o n a l d e p e n d e n c e ; the t r a n s p o r t p r o f i l e is skewed f r o m the a i r w a y center line in d i f f e r e n t ways with e a c h direct i o n o f flow t h r o u g h a bifurcation. Thus, the p e n e t r a t i o n o f gas a n d a e r o s o l particles to t h e p e r i p h e r a l airways o f the lung b y n o r m a l v e n t i l a t i o n does i n v o l v e b i f u r c a t i n g f l o w - f i e l d g e o m e t r y ; to c o m p l e t e l y u n d e r s t a n d this p r o c e s s , the simple t a p e r e d t u b e m o d e l is p r o b a b l y n o t a d e q u a t e . E v e n t u a l l y , T E M P E S T s h o u l d be a b l e to c o m p u t e t r a n s p o r t in b i f u r c a t i n g flow fields t h a t b e t t e r m o d e l a i r w a y s o f t h e lungs.
REFERENCES 1. Briant, J.K.; Lippmann, M. Particle transport through a hollow canine airway cast by high-frequency oscillatory ventilation. Exp. Lung Res. 18:385-407; 1992. 2. Gaver III, D.P.; Grotberg, J.B. An experimental investigation of oscillating flow in a tapered channel. J. Fluid Mech. 172:47-61; 1986. 3. Godleski, D.A.; Grotberg, J.B. Convection-diffusion interaction for oscillatory flow in a tapered tube. J. Biomech. Eng. 110:283-291; 1988. 4. Grotberg, J.B. Volume-cycled oscillatory flow in a tapered channel. J. Fluid Mech. 141:249-264; 1984. 5. Happel, J.; Brenner, H. Low Reynolds number hydrodynamics, with special applications to particulate media. The Hague, The Netherlands: Martinus Nijhoff Publishers; 1983. 6. Horsfield, K.; Cumming, G. Angles of branching and diameters of branches in the human bronchial tree. Bull. Math. Biophys. 29:245-259; 1967. 7. Schreck, R.M. Laminar flow through bifurcations with applications to the human lung. Ph.D. Dissertation. Evanston, IL: Northwestern University; 1972. 8. Weibel, E.R. Morphometry of the human lung. New York: Academic Press; 1963.
