ARTICLE IN PRESS
Journal of Biomechanics 37 (2004) 661–667
Modelling of peak-flow wall shear stress in major airways of the lung A.S. Green* Department of Mechanical and Design Engineering, University of Portsmouth, Anglesea Building, Anglesea Road, Portsmouth PO1 3DJ, UK Accepted 18 September 2003
Abstract Some respiratory diseases result in the inflammation of the lung airway epithelium. An associated chronic cough, as found in many cases of asthma and in long-term smokers, can exacerbate damage to the epithelial layer. It has been proposed that wall shear stresses, created by peak expiratory flow-rates during a coughing episode, are responsible. The work here uses a computational fluid dynamics technique to model peak expiratory flow in the trachea and major lung bronchi. Calculated wall shear stress values are compared to a limited set of published measurements taken from a physical model. The measurements are discussed in the context of a flow study of a complex bronchial network. A more complete picture is achieved by the calculation method, indicating, in some cases, higher maximum wall shear stresses than measured, confirming the original findings of the experimental work. Recommendations are made as to where further work would be beneficial to medical applications. r 2003 Elsevier Ltd. All rights reserved. Keywords: Wall shear stress; Lung airways; Expiration; Cough; CFD modelling
1. Introduction Asthma combines airway obstruction, increased responsiveness and inflammation, Djukanovic! et al. (1990). Clinical studies have found that chronic cough, in association with mild cases of asthma in children, may lead to exacerbation of the condition, Chang et al. (2002). Other studies show that children, predisposed to lower respiratory tract illnesses and chronic cough, may later develop epithelial damage and airway hyperresponsiveness, Heino et al. (1990). The exact effect of cough on asthma is still unknown, Young et al. (1991). More generally, adults with chronic bronchitis and productive cough, whether or not through long-term smoking, also develop inflammation of the bronchial epithelium, Cotran et al. (1999). It has been proposed that wall shear stresses, created by the peak expiratory flow-rates during cough, are responsible for additional trauma to an already inflamed epithelial layer (Cotran et al., 1999; Chowdhary et al., 1999). However, the detailed effects of shear stress, including the level at which bronchi wall epithelium sustains damage, are currently unknown. *Tel.: +44-2392-842375; fax: +44-2392-842351. E-mail address:
[email protected] (A.S. Green). 0021-9290/$ - see front matter r 2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.jbiomech.2003.09.024
In the absence of in vivo data the authors (Chowdhary et al., 1999), proposing the cough damage mechanism, undertook an experimental investigation. Wall shear stresses were measured in a physical model of the human trachea and major bronchi to assess likely magnitudes. The physical model was comprised of straight, circular cross-section, glass tubes, fused together, in an arrangement based on a well-known and adopted specification for the first three bronchus generations, Horsfield et al. (1971), shown in Table 1. The arrangement consisted of the trachea, the left and right main bronchus, the left upper and lower bronchi and the right upper, middle and lower bronchi lying in a single two-dimensional plane (Fig. 1). Note that ‘‘left’’ and ‘‘right’’ refer to the anatomical orientation of the human body. In the experiments a steady, metered air flow-rate to each of the five bronchi was established and the wall shear stresses were measured at six locations The range of airflow-rates, exhausting through the trachea was from 1 l/s, a normal tidal expiratory flow, up to a peak flow of 8 l/s representing cough. These flow-rates produce Reynolds numbers of 5394 and 43,149, respectively. Typical air velocities in the trachea, for cough, are in the Reynolds number range 104–105, Ross et al. (1955). The wall shear stresses were measured, over the range of flow-rates, using Preston tubes, Preston (1954).
ARTICLE IN PRESS A.S. Green / Journal of Biomechanics 37 (2004) 661–667
662
Nomenclature cf Rer um
skin friction coefficient, tw =ð12ru2r Þ recirculation Reynolds number, ðrur yn Þ=m mean velocity (m/s)
Table 1 Dimensions of a three-generation bronchus network from [6] Network bronchus
Length (mm)
Diameter (mm)
Angle ( )
Flow (%)
Trachea Left main Right main Left upper Left lower Right secondary Right upper Right lower Right middle
100 50 22 16 11 26 16 8 21
16 12 11 7 8 9 7 6 5
0 73 35 48 44 15 63 15 61
100 45 55 20 25 35 20 20 15
Fig. 1. The physical model of the first three bronchus generations showing the measuring locations and flow directions.
The computational fluid dynamics (CFD) technique has been used previously to study flows in a single airway bifurcation, i.e. one parent branch and two daughter branches, as well as multiple junctions in twodimensions. Most of this work has been for a system of central or smaller airway junctions within the lung and subject to laminar flow. A more realistic, three-dimensional CFD multiple junction model, in laminar flow, by Liu et al. (2002) has been validated with the velocity profile measurements of Zhao and Lieber (1994).
ur yn m r tw
maximum recirculation velocity (m/s) distance from wall to calculation node (m) dynamic viscosity (Pa s) density (kg/m3) wall shear stress (Pa)
The flow within the trachea and major bronchi is turbulent for the normal range of flow-rates. Some recent CFD work on a multiple junction model, employing more realistic, curved bronchi, compared transient velocity profiles with a validated single bifurcation, Calay et al. (2002). In this work the carina for each junction, the wall section defining the region between each daughter branch, was assumed to be sharply defined. Although this gives an unambiguous representation for modelling purposes, in reality the carina is not sharp but has a small but significant radius. The actual shape is difficult to model accurately in the absence of precise data; shear stress calculations were not reported. Early experimental studies (Chang and El Masry, 1982; Isabey and Chang, 1982), based on a 3:1, curved bronchi, physical model, presented velocity profiles for both inspiratory and expiratory flows. The expiratory flow-rates used were 0.4 and 1.7 l/s with Reynolds numbers, based on the trachea diameter, of 2123 and 8846, respectively. Shear stress measurements were not reported. The work described here uses CFD to calculate the flow structure, and particularly the wall shear stresses, within a multiple junction model for turbulent flow conditions. A geometrical representation of the experimental physical model (Chowdhary et al., 1999), described above, is used. In the absence of measured velocity profile data, the calculated mean velocity and turbulence intensity profiles are discussed in the context of other (very few) published measurements. The calculated wall shear stresses are compared to those measured in the experiments (Chowdhary et al., 1999). The implications of using straight tubes to represent the bronchi and the consequent effect on the airflow patterns are also discussed. Recommendations are made where further work would be beneficial to medical applications.
2. The CFD model CFD is a technique used to simulate the flow within or around an object of interest. The physical aspects of any fluid flow are governed by the fundamental principles of mass, momentum and energy conservation. These fundamental principles can be expressed, in their most general form, in terms of non-linear partial
ARTICLE IN PRESS A.S. Green / Journal of Biomechanics 37 (2004) 661–667
values shown in Table 1. In the calculated work it was assumed that flow entered each bronchus normally to the inlet plane, without secondary flow. At the trachea ‘‘outlet’’ a zero-gradient boundary condition was imposed.
3. Results Neither flow visualisation nor velocity profiles were reported for the experimental work (Chowdhary et al., 1999). Although uncorroborated, it is useful to gain some insight into the general flow pattern by studying the calculated velocity profiles at various locations. 3.1. Results for the 1 l/s flow-rate case The mean velocity profiles for the 1 l/s flow-rate (normal tidal expiratory) case and the full-scale bronchi geometry are shown in Figs. 2 and 3. Each velocity has been made non-dimensional with the outflow bulk velocity and its location non-dimensional with the local bronchus diameter. The profiles are located at the same
1.5 1.3 1.1
Flow velocity
0.9 0.7 0.5 A B C
0.3 0.1 -0.1 0
0.2
0.4
0.6
0.8
1
-0.3 -0.5 Bronchus diameter
Fig. 2. The calculated mean velocity profiles for locations A–C.
2
1.5
Flow velocity
differential equations. The CFD technique approximates and numerically solves the fluid flow equations by discretising them over the domain of interest, using a finite-volume mesh. The results of the iterative solution procedure may be conveniently manipulated and displayed graphically for analysis. The commercially available CFD software package FLUENT 5.5 was used in this study. Some aspects of a flow, particularly turbulence within a complex geometry, currently have to be modelled and require some experience in their use for accurate results; the implications of this are described in detail by Peric! and Ferziger (1999). In carrying out a CFD calculation it is necessary to follow a number of steps. These include, firstly, determining the important features of the flow to be modelled, creating the model geometry and its finitevolume mesh of grid cells. The fluid properties and the boundary conditions require specification before the set of equations may be solved, as would be the case for any calculation. A solid modelling tool, GAMBIT 1.3, was used to construct the geometry and finite-volume mesh, for the trachea and bronchi model, prior to exporting it to the FLUENT solver. A discussion on meshing techniques is provided by Thompson et al. (1998). The solver uses a second-order accurate discretising scheme. The standard (high Reynolds number) k2e formulation has been the most applied engineering turbulence model for many years and has advantages of simplicity and economy (Peric! and Ferziger, 1999). The RNG model is an extension to the standard k2e model to improve application to flows with separation regions and streamline curvature. In addition, the RNG model is more able in representing low Reynolds number and near-wall flows and has, therefore, been used in this work. The geometry for the physical model was created from straight tubes joining as smoothly as possible. As regards the calculations, the objective of the work was to replicate the physical model used in the experiment, rather than an idealized human bronchus. Due to the physical model being constructed of glass tubes, fused together, a sharp carina could not be assumed. In discussion with the corresponding author (Chowdhary, 2002) as to the nature of the joining technique it was decided to adopt the ‘‘physiologically realistic model’’ used in a similar CFD study (Bala! sha! zy and Hofmann, 2001). This uses smooth transition zones between parent and daughter branches with a curved carinal ridge. The three-dimensional, finite-volume mesh uses 74,490 cells, mainly of hexahedral shape in the bronchi and tetrahedral shape in the transition regions between parent and daughter bronchi. The steady-state flow was assumed to be incompressible for each of the flow-rate conditions. Due to the flow being expiratory, each bronchus ‘‘inlet’’ required a specification of boundary conditions. In the previous, published experimental work, the flow through each bronchus was adjusted to
663
1
D E F
0.5
0 0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
-0.5 Bronchus diameter
Fig. 3. The calculated mean velocity profiles for locations D–F.
1
ARTICLE IN PRESS A.S. Green / Journal of Biomechanics 37 (2004) 661–667
The shear stress profiles, calculated around the bronchus wall circumference at the six locations, are shown in Figs. 6–8. The circumferential direction is clockwise looking from the bronchus towards the flow outlet at the end of the trachea. The calculated values around the circumference are compared to the single measured value. The error bar on the measured values
1.6 1.4
Turbulence intensity
six locations where the shear stress measurements were made in the experiment, illustrated in Fig. 1. The Reynolds number, referred to the trachea diameter, is 5394. Fig. 2 shows the mean velocity profiles for locations A, B and C in the frontal plane of the geometrical arrangement. Location A is within the trachea just downstream of the joining left main bronchus (LMB) and right main bronchus (RMB). The velocity profile indicates a region of recirculating flow on the trachea wall nearest the LMB. This is due to the abrupt nature of the joining flow from the LMB via a straight tube. This has serious implications for the realism of the flow profile at this location and for the measurement of the wall shear stress; this point is discussed below. Velocity profiles within the LMB are also shown for locations B and C; these do not exhibit any recirculation. Profile C is roughly central and quite peaked in appearance. The perpendicular profile at C (not shown) is much flatter in appearance. Profile B is rather flatter and biased towards the RMB side of the bronchus. Fig. 3 shows the mean velocity profiles for locations D and E within the RMB and F within the right secondary bronchus (RBS). All three profiles, especially F, show a bias towards the LMB due to the predominant inflow direction. A recirculation region was calculated for the flows entering from the right middle and right upper bronchi. Fig. 4 shows velocity profiles for a location, A1, three tracheal diameters downstream of location A. The profiles are in the frontal plane and perpendicular to that plane. The frontal plane profile is peaked in shape, about the centreline and biased towards the RMB side of the trachea. The perpendicular velocity profile is roughly symmetrical and much flatter, with a slight central depression. The turbulence intensity profiles are shown for locations A, B and C in Fig. 5. These illustrate maximum values where the velocity gradients are greatest, e.g. in the separation region (A), and generally near walls.
B C
1 0.8 0.6 0.4
0 0
0.2
0.4
0.6
0.8
1
Bronchus diameter
Fig. 5. The calculated turbulence intensity profiles for locations A–C.
0.35 0.3 0.25 0.2 A calc F calc A F
0.15 0.1 0.05 0 0
90
180
270
360
Bronchus circumference
Fig. 6. A comparison of calculated wall shear stress profiles and single point measurements for locations A and F.
0.35
1.4
0.3
Wall shear stress (Pa s)
1.2
Flow velocity
A
1.2
0.2
Wall shear stress (Pa s)
664
1 0.8 0.6 A1 frontal
0.4 A1 perpendicular
0.2
0.25 0.2 0.15 0.1
B calc C calc B C
0.05
0
0
0
0.2
0.4
0.6
0.8
1
Bronchus diameter
Fig. 4. The calculated mean velocity profiles at A1 in the frontal and perpendicular planes.
0
90
180
270
360
Bronchus circumference
Fig. 7. A comparison of calculated wall shear stress profiles and single point measurements for locations B and C.
ARTICLE IN PRESS A.S. Green / Journal of Biomechanics 37 (2004) 661–667 0.025
0.5
Wall shear stress (non-dimensional)
0.45
Wall shear stress (Pa s)
665
0.4 0.35 0.3 0.25 0.2 0.15
D calc E calc D E
0.1 0.05 0
0.02
A (1 l/s) A (8 l/s)
0.015
0.01
0.005
0
0
90
180
270
360
0
90
180
270
360
Bronchus circumference
Bronchus circumference
Fig. 8. A comparison of calculated wall shear stress profiles and single point measurements for locations D and E.
Fig. 9. A comparison of calculated wall shear stress profiles at location A, for tracheal flow-rates of 1 and 8 l/s.
equates to 720%, the error quoted by the authors (Chowdhary et al., 1999). The profiles for locations A and F, shown in Fig. 6, indicate that there are substantial circumferential variations. The recirculation region at A occurs between 170 and 190 and the shear stress, here, is calculated to be between 0.05 and 0.11 Pa rising sharply at 180 . These may be compared to the measured value of 0.01 Pa. The shear stress variation, for the circumferential location at F, is less but does show a substantial decrease at 0/360 , where the velocities are shown to be very low (Fig. 3). The measured value for location F, at 180 , is 0.2 Pa and may be compared to the calculated value of 0.25 Pa. Fig. 7 shows the circumferential wall shear stress profiles for the LMB, at B and C. These show marked variation with small peaks and valleys. Substantial reductions for profile C, at 0/360 and at 180 , correspond with the centrally, peaked profile for the velocity in this (frontal) plane (Fig. 2). The measured value at B (0/360 ) is 0.25 Pa and corresponds to a calculated value of only 0.01–0.02 Pa. The measured value at C (180 ) is 0.16 Pa corresponding to a calculated value of 0.11 Pa. The circumferential wall shear stress profiles for the RMB, at D and E, are shown in Fig. 8. These generally exhibit less variation than those in the LMB reducing substantially only at 0/360 due to the flow bias towards the LMB. The measured value at D is 0.37 Pa corresponding to the calculated value of 0.1–0.14 Pa. At location E (180 ), the measured and calculated values of wall shear stress are 0.4 and 0.06 Pa, respectively.
Table 2 Comparison of measured and calculated wall shear stresses at 8 l/s
3.2. Results for the 8 l/s flow-rate case The Reynolds number for the 8 l/s (peak expiratory flow for cough) flow-rate case, referred to the trachea diameter, is 43,149. The general flow pattern and profile shapes at the six locations are qualitatively similar to the 1 l/s case, described above. Recirculation regions were
Location
Measured value (Pa s)
Circumferential location ( )
Calculated value (Pa s)
Calculated range (Pa s)
A B C D E F
1.75 6.03 3.10 7.63 8.89 3.54
180 0 180 0 180 180
2.70 1.95 1.15 3.17 9.05 9.30
1.19–10.32 1.95–5.11 1.67–5.38 3.17–11.48 1.70–10.5 1.38–11.41
again calculated for the flows at location A and entering the right middle bronchus. The distributions of turbulence intensity are also qualitatively similar. The shear stress profile shapes for both flow-rate cases are also similar. Fig. 9 shows the circumferential profiles for location A. The plots are shown, for comparison, as non-dimensional, i.e. tw =ð12ru2m Þ: The actual shear stress values for the 8 l/s case, at each measuring location, are summarised in Table 2. The measured, the calculated value (at the measuring point) and the range of calculated values around the circumference are compared. The authors of the experimental work (Chowdhary et al., 1999) quote a measurement error for the measurements of 72%.
4. Discussion 4.1. General model features and flow patterns An arrangement of straight tubes undoubtedly produces an inferior physical model of the first three airway generations. Acute angled junctions set up local recirculation regions not normally experienced for a steady-state flow. The curved bronchus efficiently distributes the flow when dividing (inspiration) and combining (expiration) at bronchi junctions. It also sets
ARTICLE IN PRESS 666
A.S. Green / Journal of Biomechanics 37 (2004) 661–667
up a biased secondary flow pattern, as has been demonstrated in previous work. However, for the purposes of the work introduced in this paper, the straight-tube arrangement was the only study reporting wall shear stress measurements. The use of a quasisteady flow assumption has been discussed and justified by Scherer and Burtz (1978) by comparing transient and convection terms in the Navier–Stokes equations. 4.2. Velocity profiles Little comparison may be made to previous work as regards the velocity profiles due to the use of straight rather than curved tubes. However, a qualitative comparison may be made with the general findings of Chang and El Masry (1982). Here, one noted feature of their higher Reynolds number (8864) work is the difference between the frontal plane profile and the perpendicular plane profile. The frontal plane profile is peaked on the centreline but the perpendicular plane profile is flat, this is also shown for this work in Fig. 9. The implication is that for both curved and straight bronchi, the joining of the flows from each daughter bronchus produces a qualitatively similar flow pattern. The flow pattern implies greater mass transfer off the frontal plane. The main qualitative difference between the flow profiles in straight and curved bronchi occurs at the junction of the LMB and RMB with the trachea, location A. Here, a large region of recirculating flow is shown (Fig. 2) on the LMB side of the entrance, unlike that measured in Chang and El Masry (1982). 4.3. Wall shear stress profiles The calculated wall shear stress profiles have been compared to those measured in (Chowdhary et al., 1999). However, the comparison is not straightforward due to the use of the Preston tube measuring technique used for this particular application. In the experimental paper no flow visualisation or velocity measurements were reported and hence no knowledge of the local flow pattern at each measuring location was gained. Calibration errors by the original author (Preston, 1954), and doubts about using Preston tubes in other than neutral pressure gradients, were presented by Patel (1965). Others (Kiske et al., 1981; Gasser et al., 1993) have compared the Preston tube with more elaborate devices, e.g. sub-layer fence and wall pulsed wire, for recirculating and strong adverse pressure gradients. They conclude that applying the Preston tube well away from flow disturbances leads to errors of 72%. Within flow disturbances this error is much greater although the error depends on the application. For instance, the measured wall shear stress value of 0.01 Pa, at location A, appears to be particularly low. One reason for this
may be that the Preston tube was located in a recirculation region, as predicted by the calculations. Additionally, in the lung model experimental paper (Chowdhary et al., 1999), the authors do not report using the updated calibration of Patel (1965). Thus the error value of 720% quoted for the 1 l/s experiment is probably correct, however, the 72% error quoted for the 8 l/s appears to be too low. The turbulence modelling used in the calculations relies on the so-called universal empirical constants, based on the logarithmic ‘‘law-of-the-wall’’. Wall functions have a significant influence on, for example, wall shear stress and should not be applied to recirculation regions without some consideration of the results. This cannot be achieved particularly rigorously but measurements within similar flows may be used as a rough check on the consistency of the CFD calculation. Measurements of wall shear stress have been made in recirculation regions and the skin friction coefficient cf evaluated. It has been shown that, within such a region, cf is inversely proportional to Rer and should be in the range 0:05ðRer B50Þ20:01ðRer B1100Þ; Adams and Johnson (1988). The calculated value of skin friction coefficient at location A is 0.02(Rer B137) and 0.01(Rer B1068), for the 1 and 8 l/s case, respectively. Therefore the calculated recirculation regions are at least in accordance with empirical knowledge. Some of the calculated profiles of wall shear stress exhibit characteristic ‘‘peaks and valleys’’. This is realistic behaviour and has been noted in other experimental work using a range of measuring devices, located around the circumference of a tube (Gasser et al., 1993). Generally, the calculated values follow the trends of the measured values when related to each other. However, the quantitative agreement at a particular location is of variable quality. Some of the largest and steepest ‘‘valleys’’ in the calculated profiles occur at the circumferential location of the measuring point. In the absence of measurements (Chowdhary et al., 1999) around the circumference of the measuring location, within the model bronchi, and any idea of the flow structure, it is necessary to be careful about drawing conclusions. For example the steep valleys might be skewed in the experiment. The calculated maximum values of wall shear stress, throughout, are located in the RM bronchus, just upstream of the RM/RL junction, for both flow-rate cases. At 1 l/s the maximum value is 0.9 Pa compared to the highest measured value of 0.4 Pa. At 8 l/s the maximum value is 19 Pa compared to the highest measured value of 8.9 Pa. Thus the original proposal of the authors, that substantial shear stress values may be achieved by expiratory flow-rates equivalent to cough, is confirmed. The actual potential for epithelial damage would be the subject of additional physiological studies.
ARTICLE IN PRESS A.S. Green / Journal of Biomechanics 37 (2004) 661–667
5. Conclusions Calculations of wall shear stress in a model of the first three generations of lung bronchi have been carried out to compare with the only set of measurements yet published. The use of Preston tubes for wall shear stress measurement within the bronchial network, as used in the experimental work (Chowdhary et al., 1999), may lead to errors greater than that quoted. This and the single measurement point at each location make comparison with calculations more difficult to interpret. Generally, the calculated values follow the trends of the measured values when related to each other. However, the quantitative agreement at any particular location is of variable quality. The calculated maximum values of wall shear stress, throughout, are located just upstream of the right upper bronchus junction for both flow-rate cases. At 1 l/s the maximum value is 0.9 Pa, at 8 l/s the maximum value is 19 Pa. The original proposal of the authors performing the experimental work (Chowdhary et al., 1999), that substantial shear stress values may be achieved by expiratory flow-rates equivalent to cough, is confirmed by the calculations.
6. Recommendations Wall shear stress effects on the walls of the upper # in the inflammatory airways may have an important role process of respiratory disease. The work here has shown that there is merit in improving the experimental model and wall shear stress measuring technique. It is suggested that further experimental work should use curved pipes and a careful carina geometrical definition, or range of definitions. The use of a larger scale model with Reynolds and Dean number similarity would aid in the use of better measuring techniques, e.g. pulsed wall probes, and their application over a wider area of the bronchial network.
Acknowledgements The author is grateful to Mrs. Zo.e Edwards for assistance with the figures.
References Adams, E.A., Johnson, J.P., 1988. Flow structure in the near-wall zone of a turbulent separated flow. AIAA Journal 26 (8), 932–939.
667
Bal!ash!azy, J., Hofmann, W., 2001. Fluid dynamics and related particle deposition patterns in human airway bifurcations. In: Martonen, T. (Ed.), Medical Applications of Computer Modelling. The Respiratory System, WIT Press, 2001 (Chapter 5). Calay, R.K., Kurujareon, J., Hold^, A.E., 2002. Numerical simulation of respiratory flow patterns within human lung. Research in Physiology and Neurology 130, 201–221. Chang, H.K., El Masry, O.A., 1982. A model study of flow dynamics in human central airways. Part I: axial velocity profiles. Research in Physiology 49, 75–95. Chang, A.B., Harrhy, V.A., Simpson, J., Masters, I.B., Gibson, P.G., 2002. Cough airway inflammation and mild asthma exacerbation. Archives of Disease in Childhood 86 (4), 270–275. Chowdhary, R., Singh, V., Tattersfield, A.E., Sharma, S.D., Subir, K., Gupta, A.B., 1999. Relationship of flow and cross-sectional area to frictional stress in airway models of asthma. Journal of Asthma 36(5), 419–426. Cotran, R.S., Collins, T., Kumar, V., Schmitt, W., 1999. Robbins Pathologic Basis of Disease 6th Edition. Saunders, London. Djukanovi!c, R., Roche, W.R., Wilson, J.W., Beasley, C.R.W., Twentyman, O.P., Howarth, P.H., Holgate, S.T., 1990. Mucosal inflammation in asthma. American Review of Respiratory Disease 142, 434–457. Gasser, D., Thomann, H., Dengel, P., 1993. Comparison of four methods to measure wall shear stress in a turbulent boundary layer with separation. Experiments in Fluids 15, 27–32. Heino, M., Juntunen-Backman, K., Leijala, M., Rapola, J., Laitinen, L.A., 1990. Bronchial epithelial inflammation in children with chronic cough after early lower respiratory tract illness. Respiratory Disease 141, 428–432. Horsfield, K., Dart, G., Olsen, D.E., Filley, G.F., Cumming, G., 1971. Models of the human bronchial tree. Journal of Applied Physiology 31, 207–217. Isabey, D., Chang, H.K., 1982. A model study of flow dynamics in human central airways. Part II: secondary flow velocities. Research in Physiology 49, 97–113. Kiske, S., Vasanta Ram, V., Pfarr, K., 1981. The effect of disturbed turbulence structure on the preston-tube method of measuring wall shear stress. Aeronautical Quarterly 32, 354–367. Liu, Y., So, R.M.C., Chang, C.H., 2002. Modelling the bifurcating flow in a human lung airway. Journal of Biomechanics 35, 465–473. Patel, V.C., 1965. Calibration of preston tube and limitations on its use in pressure gradients. Journal of Fluid Mechanics 23, 185–208. Peri!c, M., Ferziger, J.H., 1999. Computational Methods for Fluid Dynamics. Springer, Berlin. Preston, J.H., 1954. The determination of turbulent skin friction by means of pitot tubes. Journal of Royal Aeronautical Society 58, 109. Ross, B.B., Gramiak, R., Rahn, H., 1955. Physical dynamics of the cough mechanism. Journal of Applied Physiology 8, 264–269. Scherer, P.W., Burtz, L., 1978. Fluid mechanical experiments relevant to coughing. J. Biomechanics 11, 183–187. Thompson, J.F., Soni, B., Weathermill, N.P., 1998. Handbook of Grid Generation. CRC, Boca Raton. Young, S., Bitsakou, H., Cari!c, D., McHardy, G.J.R., 1991. Coughing can relieve or exacerbate symptoms in asthmatic patients. Respiratory Medicine 85 (Supplement A), 7–12. Zhao, Y., Lieber, B.B., 1994. Steady inspiratory flow in a model symmetric bifurcation. ASME Journal of Biomechanical Engineering 116, 488–496.