Original Article
Wall shear stress distributions in a model of normal and constricted small airways
Proc IMechE Part H: J Engineering in Medicine 2014, Vol. 228(4) 362–370 Ó IMechE 2014 Reprints and permissions: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/0954411914527586 pih.sagepub.com
David J Evans1, Anthony S Green2 and Nicholas K Thomas2
Abstract Previous studies have highlighted flow shear stress as a possible damage mechanism for small airways, in particular those liable to constriction through disease or injury due to mechanical ventilation. Flow experiments in vitro have implicated shear stress as a relevant factor for mechanotransduction pathways with respect to airway epithelial cell function. Using computational fluid dynamics analysis, this study reports velocity profiles and calculations for wall shear stress distributions in a three-generation, asymmetric section of the small airways subjected to a steady, inspiratory flow. The results show distal variation of wall shear stress distributions due to velocity gradients on the carina side of each daughter airway branch. The maximum wall shear stresses in both normal and constricted small airways are shown to exceed those calculated using data from previous simpler one-dimensional experimental analyses. These findings have implications for lung cell flow experiments involving shear stress in the consideration of both normal airway function and pathology due to mechanotransduction mechanisms.
Keywords Small airways, bronchial epithelium, velocity profiles, shear stress modelling, mechanotransduction experiments
Date received: 27 March 2013; accepted: 7 February 2014
Introduction The small airways are defined as those airways of less than approximately 2 mm internal diameter. The dichotomous branching pattern of the tracheobronchial tree results in increasingly large numbers of airways in peripheral generations. Therefore, small airways constitute the majority of lung airways, although collectively they contribute little resistance to flow. In fact, the resistance arising from these airways is only 10%–20% of the total airway resistance, and this has led to the label of being the ‘silent zone’ of the lung.1 Accordingly, pathology can arise within the small airways prior to the manifestation of symptoms, and they are recognised as relevant to the pathophysiology of diseases such as asthma, chronic obstructive pulmonary disease (COPD) and bronchiectasis. Recent asthma research has shown that bronchoconstriction is damaging to the airway in the absence of allergic inflammatory events, and mechanical stress may be implicated in airway remodelling.2 Wall shear stress, that is, a mechanical stress, is defined as ‘the frictional force per unit surface area in the direction of flow exerted at the fluid–solid interface’.3 Shear stress– induced damage to large airway mucosa caused by air
flowing over narrowed airways has been proposed.4 Shear stress, for a combination of high inspiratory flows and airway narrowing, has also been proposed as a mechanism for airway injury in ventilators.5 Additionally, shear stress may detrimentally affect the large airway epithelium in turbulent flows but should also be considered significant at generations 8–10, especially at bifurcations.6 The behaviour of flow through a simple, single, symmetric bifurcation is established.7 Measurements have been made of the steady, inspiratory flow through a symmetrical bifurcation for Reynolds numbers in the range 518 \ Re \ 2089.8 The flow Reynolds number (Re) indicates the relative influence of momentum and viscosity in the character of the flow. Momentum is more influential in high Re flows, where turbulent 1
Department of Thoracic Medicine, Hemel Hempstead Hospital, Hemel Hempstead, UK 2 Formerly at Department of Mechanical Engineering, University of Portsmouth, Portsmouth, UK Corresponding author: David J Evans, Department of Thoracic Medicine, Hemel Hempstead Hospital, Hillfield Road, Hemel Hempstead, Hertfordshire HP2 4AD, UK. Email: [email protected]
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
Evans et al.
363
effects are prevalent. Conversely, low Reynolds number flows are laminar where viscous effects dominate, for example, in the human conducting airways, typically Re . 2000 in the trachea, 1500 . Re . 500 in the main bronchi and Re \ 200 in the small airways. Published computational fluid dynamics (CFD) calculations studied the three-dimensional flow at successive symmetrical airway bifurcations for generations 5– 7 (1600 . Re . 200)9 and demonstrated important redistributions of airflow velocity profiles at bifurcations. Airway resistance has also been calculated10 for airway morphology comprising generations 0–23.11 This work was based on a simple model for steadystate, one-dimensional flow,10 in straight, circular tubes including the effects of bifurcations.12 A time-dependent, one-dimensional straight tube, branching network model of all the conducting airway generations was designed to represent the heterogeneous bronchoconstriction of the smaller airways during ventilation.5 This was accomplished by representing the airways as statistical distributions of tube diameter and extent of constriction. Their work concentrated on the smallest airways, typically of nominal diameter 0.4 mm, reporting maximum wall shear stress amplification values. Computational work on straight tube models of airway generations 3–513 and the computed tomography (CT)-derived model of a single bifurcation generation 3–414 has included the effects of airway wall compliance. These studies showed that the (flow) wall shear stress is affected by the assumed Young’s modulus Y (i.e. the stress-to-strain ratio) of the airway wall. The earlier work13 with Y values of 0.1309 MPa (longitudinal) and 0.074 MPa (circumferential) showed a 66% effect on flow shear stress at Re = 1310. The later work14 assumed a value of Y of 0.001 MPa and calculated an 80% reduction in maximum flow shear stress (Re = 183, 475). The purpose of the study reported here is to use CFD calculations of velocity and wall shear stress in a three-dimensional model representing a section of the small airways, specifically, an asymmetric multiple junction network approximating to airway generations 10–12. Two flow cases are investigated for a steady, inspiratory flow: a normal case for normally functioning small airways and a constricted case where the airway walls have contracted, simulating a diseased section of small airways. The airway walls are assumed fixed and non-compliant. In the normal case, the low Re (=37) is assumed to be so low that fluid wall interactions are minimal.14 The constricted case assumes a diseased state, where the airway walls have become thickened and non-compliant. The calculations of wall shear stress are compared to the extrapolation of previously reported results based on simpler onedimensional flow.11 Quantifying the effects of shear stress on the airway epithelium is potentially informative for cell flow experiments studying the effects of mechanotransduction with respect to epithelial cell structure and
function. Wall shear stress is usually estimated using a Poiseuille formulation for developed laminar flow in long tubes. This study will provide evidence for reconsidering the assumptions made for shear stress magnitudes used in these types of experiment.
Method Due to the size of human small airways, their heterogeneity and general inaccessibility, geometric measurements are difficult to obtain.15,16 Therefore, rather than select arbitrary dimensions and bifurcating angles, it was assumed that a typical network of the human small airways could be represented by a reduced scale version of the asymmetric system of the major bronchi. In this context, work with in vivo CT imaging of human lungs has shown that the perimeter of a large airway may be used to estimate the perimeter of a small airway.17 Accordingly, the large airway morphology18 was scaled down to represent a section of small airways encompassing generations 10–12. Details of the geometry are shown in Figure 1(a) and Table 1. The asymmetric, coplanar, branching arrangement provides for a parent (airway) branch producing two daughter branches, for example, an ‘inlet branch’ produces a ‘left branch’ and a ‘right branch’. The left branch further produces a left upper branch (L1), a left lower branch (L2) and so on. The bifurcation half-angle for each daughter branch is conventionally referred to as the angle between her centre line and the extended centre line of her parent branch. With respect to the carina, that is, the ridge between the daughter branches, a smooth transition zone was created where a radius of curvature of 10% of the major daughter diameter was used in line with accepted airway models.19 CFD is a validated method for flow calculations in a wide range of flow processes. The commercially available CFD software package FLUENT (ANSYS Inc.) was used in this study and is briefly described in Appendix 1. The calculation solver uses a second-order accurate numerical meshing scheme. The bifurcating airway structure was represented as a set of subvolumes using the GAMBIT solid modelling tool (ANSYS Inc.). Details of the sub-volumes in the region of a bifurcation are shown in Figure 1(b). Grid refinement, in the range of 70,000–120,000 sub-volumes, ensured mesh independence. The three-dimensional, finite-volume mesh used in this study comprised 95,000 sub-volumes of tetrahedral shape. An incompressible, steady-state flow rate was assumed for breathing, equivalent to that in the trachea of 0.6 L/s, a value within ranges used elsewhere, for example, 0.5–2 L/s (Pedley et al.7), 0.27–2.16 L/s (Liu et al.9) and 0.157–0.828 L/s (Xia et al.14). Flow enters the vertical inlet branch and exits the five lower branches. In reality, the flow entering a section of small airway would have a velocity distribution influenced by the upstream airway branching pattern. However, in
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
364
Proc IMechE Part H: J Engineering in Medicine 228(4)
Figure 1. (a) Schematic representation of the small airway bifurcation model, generations 9–12 used in the study, and (b) detailed view of mesh showing calculation of sub-volumes at a bifurcation.
Table 1. Dimensions of modelled small airway network and flow distribution (unconstricted case). Network airway
Modelled airway generation
Length (mm)
Diameter (mm)
Bifurcation half-angle (°)
Flow (%)
Inlet Left Right L1 L2 R1 R2 R3 R4
9 10 10 11 11 11 11 12 12
11.6 5.80 2.55 1.86 1.28 1.86 3.02 2.44 0.93
1.86 1.39 1.28 0.81 0.93 0.81 1.04 0.58 0.70
0 73 35 48 44 63 15 61 15
100 45 55 20 25 20 35 15 20
this work, the inlet branch was made relatively long sufficient for the inlet flow to develop within an entry length. Here, the inlet branch length was double the required entry length17 for both normal and constricted cases. Thus, the flow was fully developed well upstream of the first bifurcation. A proposed approach for the distribution of flow rates within the model18 was adopted, that is, the airway flow rates were based on actual measurements of lobar volumes and derived airway flows for an even distribution within the lung. This distribution, prescribed at each outlet branch of the small airway, is shown in Table 1. In the model, it was considered that this would give a better representation of the chosen small airway geometry than, for example, a constant pressure boundary condition. Therefore, there was no assumption that each of the five outlets of the small airway section showed the same pressure relative to the
inlet branch. Nonetheless, the magnitudes of the flows in each branch were calculated to be within the proposed ranges for distal bronchi (e.g. 0.107%–0.029% compared to 0.109%–0.007% of the tracheal flow rate18). Accordingly, the flow distribution at the outlet branches is as specified in Table 1 with a zero-gradient velocity profile, that is, it is assumed that the velocity profile is unchanging in each branch axial direction, consistent with Liu et al.9 and Calay et al.,19 thus maintaining the integrity of the CFD calculation. Two airway states were modelled, one with normal airway diameters and the other with the airways in a constricted state. Based on a 10th-generation diameter and flow rate, the flow was laminar with Re = 37 for the normal case. The equivalent tracheal flow rate of 0.6 L/s was maintained for the constricted case simulating mechanical ventilation of a diseased lung with Re = 146. The level of constriction was chosen to be
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
Evans et al.
365
Figure 2. Comparison of normal airway non-dimensional velocity profiles in the plane of bifurcation: (a) Stations A and A (perpendicular plane) before the first bifurcation, (b) Stations B and C in branch L1 following the first bifurcation and (c) Stations D and E in branch R1 following the first bifurcation and Station F in branch R2 following the second bifurcation.
75% of baseline diameter. In clinical terms, this value would be considered severe but falls within the range of other studies on small airway constriction, for example, 0%–90%.5,20
Results Velocity profiles Velocity profiles for the inspiratory flows for Stations A–F are plotted in Figure 2. Each velocity profile is represented in a non-dimensional format with respect to the inlet bulk velocity, that is, the inlet flow rate divided by inlet flow area. The profiles are plotted as viewed from downstream such that in a left going branch (e.g. Stations B, C and F, Figure 1(a)), the diameter at 0 is the carina side; conversely, for D and E the diameter at 1 is the carina side. Constricted case velocity profiles for Station B, for both bifurcation and perpendicular planes, are shown in Figure 3. The velocity gradients close to a branch wall are important as they are the gradients used in determining the wall shear stress (see Appendix 1). The velocity profile at Station A, in the plane of the bifurcations, upstream of the first bifurcation, is shown in Figure 2(a). The flattened profile indicates the approach of the flow to the upstream bifurcation with a lobe at the peak facing towards the left going branch, that is, Station B. The velocity profile Aperp, perpendicular to the plane of the bifurcations (anterior at 0 and posterior at 1), shows a profile generally symmetrical about the branch axis. Noting the way in
which airway divide has an important influence on the pattern of velocity gradients downstream, the velocity profile at Station B (Figure 2(b)) shows a skewed profile towards the carina side of the daughter branch airway with a sharp velocity g/radient. On the opposite, non-carina side of the branch at B, the velocity gradient is more gradual. Further downstream, at Station C, the shape of the velocity distribution, although skewed towards the carina, is more even about the centre line. Consequently, the disparity between carina and non-carina wall velocity gradients is reduced as the flow develops in the same branch. The effects of bifurcation are endorsed by the velocity profile perpendicular to the bifurcation plane at Station B (Figure 3), which is shown to be symmetrical with a flattened central peak compared with the bifurcation plane velocity profile. The velocity profile at Station D in the right going branch (Figure 2(c)) is, again, skewed towards the carina side of the branch (non-dimensional diameter location, d* = 0.62), and this is where the higher wall velocity gradient, in this plane, occurs. On the noncarina side of the branch, the velocity gradient is more gradual. Downstream, at Station E, the shape of the velocity distribution is more even about the centre line (d* = 0.53). The Station F velocity profile (Figure 2(c)) is interesting because upstream at Station E a skew in the profile exists towards the upstream carina. On entering the successive bifurcation at R2, the carina is on the opposite side of the branch. Therefore, the skew in the velocity profile, in the plane of the bifurcation,
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
366
Proc IMechE Part H: J Engineering in Medicine 228(4)
Figure 3. (a) Constricted velocity profile Station B in the plane of bifurcation (Re = 146). Constricted velocity profile Station B perpendicular to the plane of bifurcation. Comparison calculated data9 Re = 200. (b) Comparison of normal and constricted velocity profiles at Station D showing higher near-wall velocity gradients for the constricted case.
switches from one side of the branch to the other and the resulting profile is relatively symmetrical (d* = 0.51). Velocity profiles for the constricted case show a similar trend to the normal airways, although the skew effects are more marked towards the carina side. In Figure 3(a), the constricted velocity profile, at Re = 146, is compared to a published, calculated profile at a slightly greater Re = 2009 showing the similarities of the bifurcation plane profiles and skewed nature of the flow. The constricted velocity profile Bcon perp in the perpendicular plane (anterior at 0 and posterior at 1) exhibits a more symmetrical profile.9 The normal and constricted velocity profiles at D are shown (Figure 3(b)), where the maxima have been made nondimensional to 1.0. Thus, it is possible to visualise that the constricted profile exhibits higher wall velocity gradients than the normal profile. In addition to the main velocity distributions in the direction of the branches, there are secondary flow velocity vectors within the flow. For the low Reynolds number flows considered in this study, these transverse and circumferential velocities were small. For example, for the right going airway R1 at D, viewed from upstream with the carina at C, an image of the velocity vectors is shown in Figure 4 and illustrates established contrarotating vortices, at A and B. These are caused by the change in flow direction from parent to daughter branch, moving the flow towards the carina at the centre and around the transverse plane. This flow behaviour is consistent with previous descriptions7 and
Figure 4. Secondary velocity vectors at Station D for the constricted diameter case, viewed from upstream with carina at C, showing two contra-rotating vortices A and B caused by the change in flow direction from parent to daughter branch.
calculated flow patterns.9 A similar pattern of secondary flow velocity vectors was found at the other airway stations. Maximum calculated transverse velocities are about 5% and 10% of the corresponding axial velocity maxima for both unconstricted and constricted cases,
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
Evans et al.
367 respectively, also comparable to previously reported results.8
Airway wall shear stress
Figure 5. Wall shear stress profiles at points around the airway circumference (carina at 180°) at Stations A–C (normal case) showing development of circumferential shear stress in L1.
Figure 6. Wall shear stress profiles at points around the airway circumference (carina at 180°) at Stations A, D and E (normal case) showing development of circumferential shear stress in R1 and R2.
The wall shear stress values are shown for the normal case with respect to points taken around the branch circumference at Stations A–C (Figure 5) and also for Stations A, D, E and F (Figure 6). Each value is normalised with its respective wall shear stress calculated from Poiseuille flow (see Appendix 1 and Table 2). The calculations indicate circumferential variations in wall shear stress corresponding to similar variations in nearwall velocity gradients (as reported in section ‘Airway wall shear stress’). The wall shear stress profile for Station A shows calculated stresses to exceed those calculated by Poiseuille by up to 25% and further indicates the local changes to the velocity profile, that is, an increase in the local near-wall perpendicular gradient, at this location. Downstream there is a general trend for the maximum CFD calculated stresses to exceed those calculated by Poiseuille by 50%–100%. The wall shear stress profile for Station B is over 70%–90% greater in the regions 0°–90° and 270°–360° (i.e. the carina side for the left branch, Figure 5). Conversely, between 100° and 250°, the shear stress magnitudes are much less than would be calculated by Poiseuille reflecting the reduced velocity gradient on the non-carina side of the branch wall. The variation in shear stress, for Station C, is reduced due to the development of the flow in the first daughter branch. At Station C, the levels of shear stress are calculated by CFD to be up to 25% in excess of a Poiseuille-type calculation. The carina side for the right going branch (Figure 6) is in the region 90°–270°, and here, the maximum wall shear stress values for D increase by between 50% and 100% (compared to Poiseuille). On the non-carina side of the right going branch, at Station D, in the regions 0°–90° and 270°–360°, the shear stress range is between about that calculated by Poiseuille up to an excess of 50%. Just downstream at Station E, the tendency is for the circumferential variation to reduce, as at C, but remains in excess of a Poiseuille-type calculation by
Table 2. Calculated airway wall shear stresses for each station within the section of small airway. Station Unconstricted (Pa) Minimum 3 1022 Maximum 3 1022 Poiseuille calculated values103 1022 Airway constricted by 75% (Pa) Minimum Maximum Poiseuille calculated values10
A
B
C
D
E
F
3.8 4.5 3.3
2.7 6.7 3.6
3.5 5.1 3.6
6.0 11.8 5.6
6.4 8.4 5.6
6.6 8.8 6.7
2.4 2.9 2.1
0.43 7.2 2.3
2.3 4.0 2.3
3.6 11.1 3.6
3.9 8.5 3.6
3.5 6.7 4.3
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
368
Proc IMechE Part H: J Engineering in Medicine 228(4)
Figure 7. Wall shear stress profiles at points around the airway circumference (carina at 180°) at Stations B and C (normal and constricted cases).
20%–50%. The shear stresses in the third airway generation, at Station F, indicate a maximum value on the carina side of the branch at 180° about 40% above Poiseuille. However, there is another high point at 45° with an excess of 35%. This is a consequence of the velocity profile becoming less skewed due to upstream influences (i.e. the successive opposing directions of bifurcations R and R2, discussed in section ‘Airway wall shear stress’) and the weakening of secondary velocities (as Reynolds number diminishes). In line with data for velocity profiles and near-wall velocity gradients, the patterns for wall shear stress variations are also accentuated in the constricted airway case. Figure 7 shows a comparison of normal and constricted cases for circumferential wall shear stress profiles (non-dimensional) at Stations B and C, reflecting the markedly higher wall shear stresses for the constricted case. At the non-carina side of the airway (0°–45°) and (315°–360°), the increase in wall shear stress is shown in the Station C profiles as the velocity profiles develop, as discussed in section ‘Velocity profiles’. In terms of the ratio of increased wall shear stress due to the constriction of 75% over the normal case, a Poiseuille analysis would indicate a ratio of 64 (see Appendix 1). The calculated ratios for peak wall shear stresses, shown in this work, suggest a ratio range of 76–107.
Discussion Calculated results The tendency for the location of maximum axial velocity to be skewed towards the carina side of the airway, downstream of bifurcations, is evident in this study. This finding is similar to that previously described,7 measured8 and calculated9,21 elsewhere. The effect depends upon a number of variables, the most important of which is the Reynolds number. At higher Reynolds numbers of 500–1600, the secondary velocities caused by flow through a bifurcation are strong
resulting in highly skewed velocity profiles21 and consequently high wall shear stress. In the smaller airways where Reynolds numbers are 200, or less, the velocity profiles are less skewed but it is observed that they still produce wall shear stresses well in excess of that calculated by a simpler Poiseuille calculation. Other variables that influence the velocity profile at a bifurcation are flow curvature, the bifurcation daughter angles, their diameters, lengths prior to further bifurcation and flow division. In reality, the human tracheobronchial tree is highly variable, and an airway of a certain diameter may arise within a range of generations.1 Accordingly, the results should be interpreted in the context of airway heterogeneity. The calculated CFD shear stresses are compared to values derived from another study10 investigating a multiple-airway model. The upper rows of Table 2 compare the normal case for each station in the network. The comparison study10 assumes Poiseuille flow in each airway branch and relied on an empirically derived model12 to account for the additional viscous energy dissipation at upstream bifurcations. A universal constant was used which the authors acknowledged should depend on branching angle and parent-to-daughter area ratio. Thus, the true effects of these variables were not taken into account by this simpler calculation. In the absence of other directly comparable study data, the work of Wiggs et al.10 was extended to obtain wall shear stress by using the equation derived for Poiseuille flow (see Appendix 1). This exercise showed that there is a trend for shear stress maxima to be under-predicted in these previous estimates compared with the work reported here. The constricted airway case comparison is shown in the lower rows of Table 2. Data from the CFD model show shear stresses calculated to be generally two orders of magnitude higher than for the normal case. When compared with values calculated using the simpler model of Wiggs et al.,10 generally higher shear stress maxima were found. These findings reflect the substantial variation in velocity profiles calculated using a three-dimensional flow model when compared to assumptions extrapolated from a simpler onedimensional approach. The authors of previous work5 concluded that the degree of wall shear stress amplification that the highest shear stresses were found in airways constricted 60%–70% of baseline diameter and that constriction could increase maximum levels by 50fold, although localised regions could reach 600-fold. The calculated wall shear stress distributions may also be compared to previous cell flow experiments, where small airway shear stress has been studied in relation to epithelial function. Shear stresses for normal breathing used in those experiments6,22,23 were 0.03–0.05 Pa, compared to the calculated values of 0.027–0.118 Pa in this article. Therefore, from the results, both simpler one-dimensional mathematical models and in vitro cell culture experiments may have underestimated the magnitude and biological effects of
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
Evans et al.
369
shear stress. These findings might suggest that an increase in the shear stress in a straight conduit of cultured epithelial cells by 50%–100% would allow observation of the true effects of small airway bifurcations. Future work on airway wall shear stress will entail sophisticated modelling techniques where the compliance of the airway wall is taken into account.13,14 However, more experimental work on accurate values for Young’s modulus of airway wall types is needed.
Medical implications Pathology can progress within the ‘silent zone’ of the lung without recognition, partly due to the fact that standard measures of lung function are poor at identifying impairment due to small airway disease.24 This is a clinically relevant problem given that individuals with small airway disease present at an advanced stage of their condition notwithstanding the additional difficulties faced by researchers trying to understand the key pathophysiological mechanisms of these disorders. The increased shear stress found in the model of airway narrowing is consistent with the findings of other investigators. Tissues are known to be mechanosensitive where physical effects influence structure and function. Accordingly, it is plausible that deleterious effects of increased stresses, exerted upon the airways of individuals with asthma or COPD may be relevant to the pathology of these conditions. It may be that the chronic effects of airway obstruction result in levels of airway wall shear stresses that incite further damage2 in addition to any other primary elements such as allergy and infection. To date, the exact mechanism that may link airway shear stress and the tissue pathology that characterises airway disorders, such as basement membrane thickening, peribronchial fibrosis and smooth muscle hypertrophy, has not been fully established. Further detailed studies will be required to confirm this possibility. In vitro models have shown the effects of physical stimuli on airway epithelial cells with respect to elements of intracellular signalling. Data show that shear stress stimulates the release of nucleotides by airway epithelial cells and that these molecules influence basic airway functions such as ciliary beat frequency and the volume of airway surface fluid (ASF) via the activation of apical membrane purinoceptors.22,25 These investigators emphasise the importance of rhythmic mechanical stress (i.e. shear stress) for effects on airway function and host defence. Shear stress calculated in this study for normal airways may suggest the possibility of increased importance for this mechanism.
Conclusion The velocity profiles and wall shear stress distributions for a steady inspiratory flow in a representative threedimensional asymmetric section of the small airways (normal and constricted) have been calculated. The
resulting wall shear stress distributions show variations in accordance with the increased velocity gradients on the carina side of each daughter branch. Effects for velocity profiles and wall shear stresses are accentuated by airway constriction. In addition, the maximum wall shear stresses in a normal and constricted small airway are shown to exceed those calculated using previous one-dimensional deterministic analyses. Several recent in vitro studies have used shear stress in mechanotransduction experiments on airway epithelial cells. Work in this article suggests that the shear stress effects of these studies may be an underestimate with respect to cellular functions and possible airway pathology. A better understanding of this aspect of small airway physiology (and pathology) may be pertinent to the treatment of established diseases and the progression of the same. Acknowledgements The authors are grateful to Mrs Zoe¨ Edwards for assistance with Figure 1(a). Declaration of conflicting interests The authors declare that there are no conflicts of interest related to the work in this article. Funding This research was funded by the Engineering and Physical Sciences Research Council (EPRSC) following grant applications by Drs Green and Evans under Research Contracts GR/T20441/01 and GR/R56143/01. References 1. Macklem PT and Mead M. Resistance of central and peripheral airways measured by a retrograde catheter. J Appl Physiol 1967; 22(3): 395–401. 2. Grainge CL, Lau LC, Dulay V, et al. Effect of bronchoconstriction on airway remodelling in asthma. N Engl J Med 2011; 364(21): 2006–2015. 3. Liu M, Tanswell AK and Post M. Mechanical forceinduced signal transduction in lung cells. Am J Physiol 1999; 277: 667–683. 4. Chowdhary R, Singh V, Tattersfield AE, et al. Relationship of flow and cross-sectional area to frictional stress in airway models of asthma. J Asthma 1999; 36(5): 419– 426. 5. Nucci G, Suki B and Lutchen K. Modeling airflowrelated shear stress during heterogeneous constriction and mechanical ventilation. J Appl Physiol 2003; 95(1): 348–356. 6. Kotani M, Kotani T, Li Z, et al. Reduced inspiratory flow attenuates IL-8 release and MAPK activation of lung overstretch. Eur Respir J 2004; 24: 238–246. 7. Pedley TJ. Pulmonary fluid dynamics. Annu Rev Fluid Mech 1977; 9: 229–274. 8. Zhao Y and Lieber BB. Steady inspiratory flow in a model symmetric bifurcation. J Biomech Eng 1994; 116(4): 488–496.
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
370
Proc IMechE Part H: J Engineering in Medicine 228(4)
9. Liu Y, So RMC and Zhang CH. Modeling the bifurcating flow in a human lung airway. J Biomech 2002; 35: 465–473. 10. Wiggs BR, Bosken C, Pare PD, et al. A model of airway narrowing in asthma and in chronic obstructive pulmonary disease. Am Rev Respir Dis 1992; 145: 1251–1258. 11. Weibel ER. Morphometry of the human lung. New York: Springer-Verlag, 1963. 12. Pedley TJ, Schroter RC and Sudlow MF. Energy losses and pressure drop in models of human airways. Respir Physiol 1970; 9: 371–386. 13. Koombua K and Pidaparti RM. Inhalation induced stresses and flow characteristics in human airways through fluid–structure interaction analysis. Model Simulat Eng 2008; 358748 (8 pp.), http://dx.doi.org/10.1155/2008/ 358748 14. Xia G, Tawhai M, Hoffman E, et al. Airway wall stiffening increases peak wall shear stress: a fluid–structure interaction study in rigid and compliant airways. Ann Biomed Eng 2010; 38(5): 1836–1853. 15. Philips CG and Kaye SR. On the asymmetry of bifurcations in the bronchial tree. Respir Physiol 1997; 107: 85–98. 16. Nakano Y, Wong JC, de Jong PA, et al. The prediction of small airway dimensions using computed tomography. Am J Respir Crit Care Med 2005; 171: 142–146. 17. Durst F, Ray S, Unsal B, et al. The development lengths of laminar pipe and channel flows. J Fluid Eng: T ASME 2005; 127: 1154–1160. 18. Horsfield K, Dart G, Olsen DE, et al. Models of the human bronchial tree. J Appl Physiol 1971; 31: 207–217. 19. Calay RK, Kurujareon J and Holdø AE. Numerical simulation of respiratory flow patterns within a human lung. Resp Physiol Neurobiol 2002; 130: 201–221. 20. Tgavalekos NT, Tawhai M, Scott Harris R, et al. Identifying airways responsible for heterogeneous ventilation and mechanical dysfunction in asthma: an image functional modeling approach. J Appl Physiol 2005; 99: 2388–2397. 21. Mauroy B, Filoche M, Andrade J, et al. Interplay between geometry and flow distribution in an airway tree. Phys Rev Lett 2003; 90(14): 148101. 22. Tarran R, Button B, Picher M, et al. Normal and cystic fibrosis airway surface liquid homeostasis – the effects of phasic shear stress and viral infections. J Biol Chem 2005; 280(42): 35751–35759. 23. Sidhaye VK, Schweitzer KS, Caterina MJ, et al. Shear stress regulates aquaporin-5 and airway epithelial barrier function. Proc Natl Acad Sci USA 2008; 105(9): 3345–3350. 24. Evans DJ and Green M. Small airways, a time to revisit? Thorax 1998; 53: 629–630. 25. Button B, Picher M and Boucher RC. Differential effects of cyclic and constant stress on ATP release and
mucociliary transport by human airway epithelia. J Physiol 2007; 580(2): 577–592.
Appendix 1 The CFD model A CFD model has been chosen to examine velocity profiles and shear stress distributions in a model of small airways. The CFD technique approximates and numerically solves the incompressible Navier–Stokes fluid flow equations (1) and (2). There is no analytical solution for these equations in a flow system due to the complexity of that geometry (such as a bifurcating network of tubes). Therefore, the equations have to be solved numerically using a computer software program, and such techniques have been developed and refined over the last 40 years r u=0 ru ru = rp + mr2 u
ð1Þ ð2Þ
where u is the velocity vector (m/s), p is the pressure (Pa), r is the air density (kg/m3) and m is the air molecular viscosity (Pa s).
Wall shear stress The wall shear stress t only occurs where flow interacts with a surface (i.e. an airway wall), and it is a derived quantity defined as the product of molecular viscosity with the local perpendicular velocity gradient, t = m(du/dy). Poiseuille flow assumes ideal conditions of a long straight pipe18 of constant circular diameter and a steady, laminar flow where the velocity profile has developed to a parabolic velocity shape; velocity is inversely proportional to the square of the pipe diameter. The wall shear stress may be calculated5,14 as t=
32 mQ pd3
where Q is the airway volumetric flow rate (m3/s) and d is the airway diameter (m), and thus, the wall shear stress is inversely proportional to the cube of the pipe diameter. The degree to which conditions depart from the ideal, as in a flow bifurcation with pipe lengths of very few diameters, may be readily compared by nondimensionalising wall shear stresses to those derived from Poiseuille.
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
Wall shear stress distributions in a model of normal and constricted small airways
Proc IMechE Part H: J Engineering in Medicine 2014, Vol. 228(4) 362–370 Ó IMechE 2014 Reprints and permissions: sagepub.co.uk/journalsPermissions.nav DOI: 10.1177/0954411914527586 pih.sagepub.com
David J Evans1, Anthony S Green2 and Nicholas K Thomas2
Abstract Previous studies have highlighted flow shear stress as a possible damage mechanism for small airways, in particular those liable to constriction through disease or injury due to mechanical ventilation. Flow experiments in vitro have implicated shear stress as a relevant factor for mechanotransduction pathways with respect to airway epithelial cell function. Using computational fluid dynamics analysis, this study reports velocity profiles and calculations for wall shear stress distributions in a three-generation, asymmetric section of the small airways subjected to a steady, inspiratory flow. The results show distal variation of wall shear stress distributions due to velocity gradients on the carina side of each daughter airway branch. The maximum wall shear stresses in both normal and constricted small airways are shown to exceed those calculated using data from previous simpler one-dimensional experimental analyses. These findings have implications for lung cell flow experiments involving shear stress in the consideration of both normal airway function and pathology due to mechanotransduction mechanisms.
Keywords Small airways, bronchial epithelium, velocity profiles, shear stress modelling, mechanotransduction experiments
Date received: 27 March 2013; accepted: 7 February 2014
Introduction The small airways are defined as those airways of less than approximately 2 mm internal diameter. The dichotomous branching pattern of the tracheobronchial tree results in increasingly large numbers of airways in peripheral generations. Therefore, small airways constitute the majority of lung airways, although collectively they contribute little resistance to flow. In fact, the resistance arising from these airways is only 10%–20% of the total airway resistance, and this has led to the label of being the ‘silent zone’ of the lung.1 Accordingly, pathology can arise within the small airways prior to the manifestation of symptoms, and they are recognised as relevant to the pathophysiology of diseases such as asthma, chronic obstructive pulmonary disease (COPD) and bronchiectasis. Recent asthma research has shown that bronchoconstriction is damaging to the airway in the absence of allergic inflammatory events, and mechanical stress may be implicated in airway remodelling.2 Wall shear stress, that is, a mechanical stress, is defined as ‘the frictional force per unit surface area in the direction of flow exerted at the fluid–solid interface’.3 Shear stress– induced damage to large airway mucosa caused by air
flowing over narrowed airways has been proposed.4 Shear stress, for a combination of high inspiratory flows and airway narrowing, has also been proposed as a mechanism for airway injury in ventilators.5 Additionally, shear stress may detrimentally affect the large airway epithelium in turbulent flows but should also be considered significant at generations 8–10, especially at bifurcations.6 The behaviour of flow through a simple, single, symmetric bifurcation is established.7 Measurements have been made of the steady, inspiratory flow through a symmetrical bifurcation for Reynolds numbers in the range 518 \ Re \ 2089.8 The flow Reynolds number (Re) indicates the relative influence of momentum and viscosity in the character of the flow. Momentum is more influential in high Re flows, where turbulent 1
Department of Thoracic Medicine, Hemel Hempstead Hospital, Hemel Hempstead, UK 2 Formerly at Department of Mechanical Engineering, University of Portsmouth, Portsmouth, UK Corresponding author: David J Evans, Department of Thoracic Medicine, Hemel Hempstead Hospital, Hillfield Road, Hemel Hempstead, Hertfordshire HP2 4AD, UK. Email: [email protected]
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
Evans et al.
363
effects are prevalent. Conversely, low Reynolds number flows are laminar where viscous effects dominate, for example, in the human conducting airways, typically Re . 2000 in the trachea, 1500 . Re . 500 in the main bronchi and Re \ 200 in the small airways. Published computational fluid dynamics (CFD) calculations studied the three-dimensional flow at successive symmetrical airway bifurcations for generations 5– 7 (1600 . Re . 200)9 and demonstrated important redistributions of airflow velocity profiles at bifurcations. Airway resistance has also been calculated10 for airway morphology comprising generations 0–23.11 This work was based on a simple model for steadystate, one-dimensional flow,10 in straight, circular tubes including the effects of bifurcations.12 A time-dependent, one-dimensional straight tube, branching network model of all the conducting airway generations was designed to represent the heterogeneous bronchoconstriction of the smaller airways during ventilation.5 This was accomplished by representing the airways as statistical distributions of tube diameter and extent of constriction. Their work concentrated on the smallest airways, typically of nominal diameter 0.4 mm, reporting maximum wall shear stress amplification values. Computational work on straight tube models of airway generations 3–513 and the computed tomography (CT)-derived model of a single bifurcation generation 3–414 has included the effects of airway wall compliance. These studies showed that the (flow) wall shear stress is affected by the assumed Young’s modulus Y (i.e. the stress-to-strain ratio) of the airway wall. The earlier work13 with Y values of 0.1309 MPa (longitudinal) and 0.074 MPa (circumferential) showed a 66% effect on flow shear stress at Re = 1310. The later work14 assumed a value of Y of 0.001 MPa and calculated an 80% reduction in maximum flow shear stress (Re = 183, 475). The purpose of the study reported here is to use CFD calculations of velocity and wall shear stress in a three-dimensional model representing a section of the small airways, specifically, an asymmetric multiple junction network approximating to airway generations 10–12. Two flow cases are investigated for a steady, inspiratory flow: a normal case for normally functioning small airways and a constricted case where the airway walls have contracted, simulating a diseased section of small airways. The airway walls are assumed fixed and non-compliant. In the normal case, the low Re (=37) is assumed to be so low that fluid wall interactions are minimal.14 The constricted case assumes a diseased state, where the airway walls have become thickened and non-compliant. The calculations of wall shear stress are compared to the extrapolation of previously reported results based on simpler onedimensional flow.11 Quantifying the effects of shear stress on the airway epithelium is potentially informative for cell flow experiments studying the effects of mechanotransduction with respect to epithelial cell structure and
function. Wall shear stress is usually estimated using a Poiseuille formulation for developed laminar flow in long tubes. This study will provide evidence for reconsidering the assumptions made for shear stress magnitudes used in these types of experiment.
Method Due to the size of human small airways, their heterogeneity and general inaccessibility, geometric measurements are difficult to obtain.15,16 Therefore, rather than select arbitrary dimensions and bifurcating angles, it was assumed that a typical network of the human small airways could be represented by a reduced scale version of the asymmetric system of the major bronchi. In this context, work with in vivo CT imaging of human lungs has shown that the perimeter of a large airway may be used to estimate the perimeter of a small airway.17 Accordingly, the large airway morphology18 was scaled down to represent a section of small airways encompassing generations 10–12. Details of the geometry are shown in Figure 1(a) and Table 1. The asymmetric, coplanar, branching arrangement provides for a parent (airway) branch producing two daughter branches, for example, an ‘inlet branch’ produces a ‘left branch’ and a ‘right branch’. The left branch further produces a left upper branch (L1), a left lower branch (L2) and so on. The bifurcation half-angle for each daughter branch is conventionally referred to as the angle between her centre line and the extended centre line of her parent branch. With respect to the carina, that is, the ridge between the daughter branches, a smooth transition zone was created where a radius of curvature of 10% of the major daughter diameter was used in line with accepted airway models.19 CFD is a validated method for flow calculations in a wide range of flow processes. The commercially available CFD software package FLUENT (ANSYS Inc.) was used in this study and is briefly described in Appendix 1. The calculation solver uses a second-order accurate numerical meshing scheme. The bifurcating airway structure was represented as a set of subvolumes using the GAMBIT solid modelling tool (ANSYS Inc.). Details of the sub-volumes in the region of a bifurcation are shown in Figure 1(b). Grid refinement, in the range of 70,000–120,000 sub-volumes, ensured mesh independence. The three-dimensional, finite-volume mesh used in this study comprised 95,000 sub-volumes of tetrahedral shape. An incompressible, steady-state flow rate was assumed for breathing, equivalent to that in the trachea of 0.6 L/s, a value within ranges used elsewhere, for example, 0.5–2 L/s (Pedley et al.7), 0.27–2.16 L/s (Liu et al.9) and 0.157–0.828 L/s (Xia et al.14). Flow enters the vertical inlet branch and exits the five lower branches. In reality, the flow entering a section of small airway would have a velocity distribution influenced by the upstream airway branching pattern. However, in
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
364
Proc IMechE Part H: J Engineering in Medicine 228(4)
Figure 1. (a) Schematic representation of the small airway bifurcation model, generations 9–12 used in the study, and (b) detailed view of mesh showing calculation of sub-volumes at a bifurcation.
Table 1. Dimensions of modelled small airway network and flow distribution (unconstricted case). Network airway
Modelled airway generation
Length (mm)
Diameter (mm)
Bifurcation half-angle (°)
Flow (%)
Inlet Left Right L1 L2 R1 R2 R3 R4
9 10 10 11 11 11 11 12 12
11.6 5.80 2.55 1.86 1.28 1.86 3.02 2.44 0.93
1.86 1.39 1.28 0.81 0.93 0.81 1.04 0.58 0.70
0 73 35 48 44 63 15 61 15
100 45 55 20 25 20 35 15 20
this work, the inlet branch was made relatively long sufficient for the inlet flow to develop within an entry length. Here, the inlet branch length was double the required entry length17 for both normal and constricted cases. Thus, the flow was fully developed well upstream of the first bifurcation. A proposed approach for the distribution of flow rates within the model18 was adopted, that is, the airway flow rates were based on actual measurements of lobar volumes and derived airway flows for an even distribution within the lung. This distribution, prescribed at each outlet branch of the small airway, is shown in Table 1. In the model, it was considered that this would give a better representation of the chosen small airway geometry than, for example, a constant pressure boundary condition. Therefore, there was no assumption that each of the five outlets of the small airway section showed the same pressure relative to the
inlet branch. Nonetheless, the magnitudes of the flows in each branch were calculated to be within the proposed ranges for distal bronchi (e.g. 0.107%–0.029% compared to 0.109%–0.007% of the tracheal flow rate18). Accordingly, the flow distribution at the outlet branches is as specified in Table 1 with a zero-gradient velocity profile, that is, it is assumed that the velocity profile is unchanging in each branch axial direction, consistent with Liu et al.9 and Calay et al.,19 thus maintaining the integrity of the CFD calculation. Two airway states were modelled, one with normal airway diameters and the other with the airways in a constricted state. Based on a 10th-generation diameter and flow rate, the flow was laminar with Re = 37 for the normal case. The equivalent tracheal flow rate of 0.6 L/s was maintained for the constricted case simulating mechanical ventilation of a diseased lung with Re = 146. The level of constriction was chosen to be
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
Evans et al.
365
Figure 2. Comparison of normal airway non-dimensional velocity profiles in the plane of bifurcation: (a) Stations A and A (perpendicular plane) before the first bifurcation, (b) Stations B and C in branch L1 following the first bifurcation and (c) Stations D and E in branch R1 following the first bifurcation and Station F in branch R2 following the second bifurcation.
75% of baseline diameter. In clinical terms, this value would be considered severe but falls within the range of other studies on small airway constriction, for example, 0%–90%.5,20
Results Velocity profiles Velocity profiles for the inspiratory flows for Stations A–F are plotted in Figure 2. Each velocity profile is represented in a non-dimensional format with respect to the inlet bulk velocity, that is, the inlet flow rate divided by inlet flow area. The profiles are plotted as viewed from downstream such that in a left going branch (e.g. Stations B, C and F, Figure 1(a)), the diameter at 0 is the carina side; conversely, for D and E the diameter at 1 is the carina side. Constricted case velocity profiles for Station B, for both bifurcation and perpendicular planes, are shown in Figure 3. The velocity gradients close to a branch wall are important as they are the gradients used in determining the wall shear stress (see Appendix 1). The velocity profile at Station A, in the plane of the bifurcations, upstream of the first bifurcation, is shown in Figure 2(a). The flattened profile indicates the approach of the flow to the upstream bifurcation with a lobe at the peak facing towards the left going branch, that is, Station B. The velocity profile Aperp, perpendicular to the plane of the bifurcations (anterior at 0 and posterior at 1), shows a profile generally symmetrical about the branch axis. Noting the way in
which airway divide has an important influence on the pattern of velocity gradients downstream, the velocity profile at Station B (Figure 2(b)) shows a skewed profile towards the carina side of the daughter branch airway with a sharp velocity g/radient. On the opposite, non-carina side of the branch at B, the velocity gradient is more gradual. Further downstream, at Station C, the shape of the velocity distribution, although skewed towards the carina, is more even about the centre line. Consequently, the disparity between carina and non-carina wall velocity gradients is reduced as the flow develops in the same branch. The effects of bifurcation are endorsed by the velocity profile perpendicular to the bifurcation plane at Station B (Figure 3), which is shown to be symmetrical with a flattened central peak compared with the bifurcation plane velocity profile. The velocity profile at Station D in the right going branch (Figure 2(c)) is, again, skewed towards the carina side of the branch (non-dimensional diameter location, d* = 0.62), and this is where the higher wall velocity gradient, in this plane, occurs. On the noncarina side of the branch, the velocity gradient is more gradual. Downstream, at Station E, the shape of the velocity distribution is more even about the centre line (d* = 0.53). The Station F velocity profile (Figure 2(c)) is interesting because upstream at Station E a skew in the profile exists towards the upstream carina. On entering the successive bifurcation at R2, the carina is on the opposite side of the branch. Therefore, the skew in the velocity profile, in the plane of the bifurcation,
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
366
Proc IMechE Part H: J Engineering in Medicine 228(4)
Figure 3. (a) Constricted velocity profile Station B in the plane of bifurcation (Re = 146). Constricted velocity profile Station B perpendicular to the plane of bifurcation. Comparison calculated data9 Re = 200. (b) Comparison of normal and constricted velocity profiles at Station D showing higher near-wall velocity gradients for the constricted case.
switches from one side of the branch to the other and the resulting profile is relatively symmetrical (d* = 0.51). Velocity profiles for the constricted case show a similar trend to the normal airways, although the skew effects are more marked towards the carina side. In Figure 3(a), the constricted velocity profile, at Re = 146, is compared to a published, calculated profile at a slightly greater Re = 2009 showing the similarities of the bifurcation plane profiles and skewed nature of the flow. The constricted velocity profile Bcon perp in the perpendicular plane (anterior at 0 and posterior at 1) exhibits a more symmetrical profile.9 The normal and constricted velocity profiles at D are shown (Figure 3(b)), where the maxima have been made nondimensional to 1.0. Thus, it is possible to visualise that the constricted profile exhibits higher wall velocity gradients than the normal profile. In addition to the main velocity distributions in the direction of the branches, there are secondary flow velocity vectors within the flow. For the low Reynolds number flows considered in this study, these transverse and circumferential velocities were small. For example, for the right going airway R1 at D, viewed from upstream with the carina at C, an image of the velocity vectors is shown in Figure 4 and illustrates established contrarotating vortices, at A and B. These are caused by the change in flow direction from parent to daughter branch, moving the flow towards the carina at the centre and around the transverse plane. This flow behaviour is consistent with previous descriptions7 and
Figure 4. Secondary velocity vectors at Station D for the constricted diameter case, viewed from upstream with carina at C, showing two contra-rotating vortices A and B caused by the change in flow direction from parent to daughter branch.
calculated flow patterns.9 A similar pattern of secondary flow velocity vectors was found at the other airway stations. Maximum calculated transverse velocities are about 5% and 10% of the corresponding axial velocity maxima for both unconstricted and constricted cases,
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
Evans et al.
367 respectively, also comparable to previously reported results.8
Airway wall shear stress
Figure 5. Wall shear stress profiles at points around the airway circumference (carina at 180°) at Stations A–C (normal case) showing development of circumferential shear stress in L1.
Figure 6. Wall shear stress profiles at points around the airway circumference (carina at 180°) at Stations A, D and E (normal case) showing development of circumferential shear stress in R1 and R2.
The wall shear stress values are shown for the normal case with respect to points taken around the branch circumference at Stations A–C (Figure 5) and also for Stations A, D, E and F (Figure 6). Each value is normalised with its respective wall shear stress calculated from Poiseuille flow (see Appendix 1 and Table 2). The calculations indicate circumferential variations in wall shear stress corresponding to similar variations in nearwall velocity gradients (as reported in section ‘Airway wall shear stress’). The wall shear stress profile for Station A shows calculated stresses to exceed those calculated by Poiseuille by up to 25% and further indicates the local changes to the velocity profile, that is, an increase in the local near-wall perpendicular gradient, at this location. Downstream there is a general trend for the maximum CFD calculated stresses to exceed those calculated by Poiseuille by 50%–100%. The wall shear stress profile for Station B is over 70%–90% greater in the regions 0°–90° and 270°–360° (i.e. the carina side for the left branch, Figure 5). Conversely, between 100° and 250°, the shear stress magnitudes are much less than would be calculated by Poiseuille reflecting the reduced velocity gradient on the non-carina side of the branch wall. The variation in shear stress, for Station C, is reduced due to the development of the flow in the first daughter branch. At Station C, the levels of shear stress are calculated by CFD to be up to 25% in excess of a Poiseuille-type calculation. The carina side for the right going branch (Figure 6) is in the region 90°–270°, and here, the maximum wall shear stress values for D increase by between 50% and 100% (compared to Poiseuille). On the non-carina side of the right going branch, at Station D, in the regions 0°–90° and 270°–360°, the shear stress range is between about that calculated by Poiseuille up to an excess of 50%. Just downstream at Station E, the tendency is for the circumferential variation to reduce, as at C, but remains in excess of a Poiseuille-type calculation by
Table 2. Calculated airway wall shear stresses for each station within the section of small airway. Station Unconstricted (Pa) Minimum 3 1022 Maximum 3 1022 Poiseuille calculated values103 1022 Airway constricted by 75% (Pa) Minimum Maximum Poiseuille calculated values10
A
B
C
D
E
F
3.8 4.5 3.3
2.7 6.7 3.6
3.5 5.1 3.6
6.0 11.8 5.6
6.4 8.4 5.6
6.6 8.8 6.7
2.4 2.9 2.1
0.43 7.2 2.3
2.3 4.0 2.3
3.6 11.1 3.6
3.9 8.5 3.6
3.5 6.7 4.3
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
368
Proc IMechE Part H: J Engineering in Medicine 228(4)
Figure 7. Wall shear stress profiles at points around the airway circumference (carina at 180°) at Stations B and C (normal and constricted cases).
20%–50%. The shear stresses in the third airway generation, at Station F, indicate a maximum value on the carina side of the branch at 180° about 40% above Poiseuille. However, there is another high point at 45° with an excess of 35%. This is a consequence of the velocity profile becoming less skewed due to upstream influences (i.e. the successive opposing directions of bifurcations R and R2, discussed in section ‘Airway wall shear stress’) and the weakening of secondary velocities (as Reynolds number diminishes). In line with data for velocity profiles and near-wall velocity gradients, the patterns for wall shear stress variations are also accentuated in the constricted airway case. Figure 7 shows a comparison of normal and constricted cases for circumferential wall shear stress profiles (non-dimensional) at Stations B and C, reflecting the markedly higher wall shear stresses for the constricted case. At the non-carina side of the airway (0°–45°) and (315°–360°), the increase in wall shear stress is shown in the Station C profiles as the velocity profiles develop, as discussed in section ‘Velocity profiles’. In terms of the ratio of increased wall shear stress due to the constriction of 75% over the normal case, a Poiseuille analysis would indicate a ratio of 64 (see Appendix 1). The calculated ratios for peak wall shear stresses, shown in this work, suggest a ratio range of 76–107.
Discussion Calculated results The tendency for the location of maximum axial velocity to be skewed towards the carina side of the airway, downstream of bifurcations, is evident in this study. This finding is similar to that previously described,7 measured8 and calculated9,21 elsewhere. The effect depends upon a number of variables, the most important of which is the Reynolds number. At higher Reynolds numbers of 500–1600, the secondary velocities caused by flow through a bifurcation are strong
resulting in highly skewed velocity profiles21 and consequently high wall shear stress. In the smaller airways where Reynolds numbers are 200, or less, the velocity profiles are less skewed but it is observed that they still produce wall shear stresses well in excess of that calculated by a simpler Poiseuille calculation. Other variables that influence the velocity profile at a bifurcation are flow curvature, the bifurcation daughter angles, their diameters, lengths prior to further bifurcation and flow division. In reality, the human tracheobronchial tree is highly variable, and an airway of a certain diameter may arise within a range of generations.1 Accordingly, the results should be interpreted in the context of airway heterogeneity. The calculated CFD shear stresses are compared to values derived from another study10 investigating a multiple-airway model. The upper rows of Table 2 compare the normal case for each station in the network. The comparison study10 assumes Poiseuille flow in each airway branch and relied on an empirically derived model12 to account for the additional viscous energy dissipation at upstream bifurcations. A universal constant was used which the authors acknowledged should depend on branching angle and parent-to-daughter area ratio. Thus, the true effects of these variables were not taken into account by this simpler calculation. In the absence of other directly comparable study data, the work of Wiggs et al.10 was extended to obtain wall shear stress by using the equation derived for Poiseuille flow (see Appendix 1). This exercise showed that there is a trend for shear stress maxima to be under-predicted in these previous estimates compared with the work reported here. The constricted airway case comparison is shown in the lower rows of Table 2. Data from the CFD model show shear stresses calculated to be generally two orders of magnitude higher than for the normal case. When compared with values calculated using the simpler model of Wiggs et al.,10 generally higher shear stress maxima were found. These findings reflect the substantial variation in velocity profiles calculated using a three-dimensional flow model when compared to assumptions extrapolated from a simpler onedimensional approach. The authors of previous work5 concluded that the degree of wall shear stress amplification that the highest shear stresses were found in airways constricted 60%–70% of baseline diameter and that constriction could increase maximum levels by 50fold, although localised regions could reach 600-fold. The calculated wall shear stress distributions may also be compared to previous cell flow experiments, where small airway shear stress has been studied in relation to epithelial function. Shear stresses for normal breathing used in those experiments6,22,23 were 0.03–0.05 Pa, compared to the calculated values of 0.027–0.118 Pa in this article. Therefore, from the results, both simpler one-dimensional mathematical models and in vitro cell culture experiments may have underestimated the magnitude and biological effects of
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
Evans et al.
369
shear stress. These findings might suggest that an increase in the shear stress in a straight conduit of cultured epithelial cells by 50%–100% would allow observation of the true effects of small airway bifurcations. Future work on airway wall shear stress will entail sophisticated modelling techniques where the compliance of the airway wall is taken into account.13,14 However, more experimental work on accurate values for Young’s modulus of airway wall types is needed.
Medical implications Pathology can progress within the ‘silent zone’ of the lung without recognition, partly due to the fact that standard measures of lung function are poor at identifying impairment due to small airway disease.24 This is a clinically relevant problem given that individuals with small airway disease present at an advanced stage of their condition notwithstanding the additional difficulties faced by researchers trying to understand the key pathophysiological mechanisms of these disorders. The increased shear stress found in the model of airway narrowing is consistent with the findings of other investigators. Tissues are known to be mechanosensitive where physical effects influence structure and function. Accordingly, it is plausible that deleterious effects of increased stresses, exerted upon the airways of individuals with asthma or COPD may be relevant to the pathology of these conditions. It may be that the chronic effects of airway obstruction result in levels of airway wall shear stresses that incite further damage2 in addition to any other primary elements such as allergy and infection. To date, the exact mechanism that may link airway shear stress and the tissue pathology that characterises airway disorders, such as basement membrane thickening, peribronchial fibrosis and smooth muscle hypertrophy, has not been fully established. Further detailed studies will be required to confirm this possibility. In vitro models have shown the effects of physical stimuli on airway epithelial cells with respect to elements of intracellular signalling. Data show that shear stress stimulates the release of nucleotides by airway epithelial cells and that these molecules influence basic airway functions such as ciliary beat frequency and the volume of airway surface fluid (ASF) via the activation of apical membrane purinoceptors.22,25 These investigators emphasise the importance of rhythmic mechanical stress (i.e. shear stress) for effects on airway function and host defence. Shear stress calculated in this study for normal airways may suggest the possibility of increased importance for this mechanism.
Conclusion The velocity profiles and wall shear stress distributions for a steady inspiratory flow in a representative threedimensional asymmetric section of the small airways (normal and constricted) have been calculated. The
resulting wall shear stress distributions show variations in accordance with the increased velocity gradients on the carina side of each daughter branch. Effects for velocity profiles and wall shear stresses are accentuated by airway constriction. In addition, the maximum wall shear stresses in a normal and constricted small airway are shown to exceed those calculated using previous one-dimensional deterministic analyses. Several recent in vitro studies have used shear stress in mechanotransduction experiments on airway epithelial cells. Work in this article suggests that the shear stress effects of these studies may be an underestimate with respect to cellular functions and possible airway pathology. A better understanding of this aspect of small airway physiology (and pathology) may be pertinent to the treatment of established diseases and the progression of the same. Acknowledgements The authors are grateful to Mrs Zoe¨ Edwards for assistance with Figure 1(a). Declaration of conflicting interests The authors declare that there are no conflicts of interest related to the work in this article. Funding This research was funded by the Engineering and Physical Sciences Research Council (EPRSC) following grant applications by Drs Green and Evans under Research Contracts GR/T20441/01 and GR/R56143/01. References 1. Macklem PT and Mead M. Resistance of central and peripheral airways measured by a retrograde catheter. J Appl Physiol 1967; 22(3): 395–401. 2. Grainge CL, Lau LC, Dulay V, et al. Effect of bronchoconstriction on airway remodelling in asthma. N Engl J Med 2011; 364(21): 2006–2015. 3. Liu M, Tanswell AK and Post M. Mechanical forceinduced signal transduction in lung cells. Am J Physiol 1999; 277: 667–683. 4. Chowdhary R, Singh V, Tattersfield AE, et al. Relationship of flow and cross-sectional area to frictional stress in airway models of asthma. J Asthma 1999; 36(5): 419– 426. 5. Nucci G, Suki B and Lutchen K. Modeling airflowrelated shear stress during heterogeneous constriction and mechanical ventilation. J Appl Physiol 2003; 95(1): 348–356. 6. Kotani M, Kotani T, Li Z, et al. Reduced inspiratory flow attenuates IL-8 release and MAPK activation of lung overstretch. Eur Respir J 2004; 24: 238–246. 7. Pedley TJ. Pulmonary fluid dynamics. Annu Rev Fluid Mech 1977; 9: 229–274. 8. Zhao Y and Lieber BB. Steady inspiratory flow in a model symmetric bifurcation. J Biomech Eng 1994; 116(4): 488–496.
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
370
Proc IMechE Part H: J Engineering in Medicine 228(4)
9. Liu Y, So RMC and Zhang CH. Modeling the bifurcating flow in a human lung airway. J Biomech 2002; 35: 465–473. 10. Wiggs BR, Bosken C, Pare PD, et al. A model of airway narrowing in asthma and in chronic obstructive pulmonary disease. Am Rev Respir Dis 1992; 145: 1251–1258. 11. Weibel ER. Morphometry of the human lung. New York: Springer-Verlag, 1963. 12. Pedley TJ, Schroter RC and Sudlow MF. Energy losses and pressure drop in models of human airways. Respir Physiol 1970; 9: 371–386. 13. Koombua K and Pidaparti RM. Inhalation induced stresses and flow characteristics in human airways through fluid–structure interaction analysis. Model Simulat Eng 2008; 358748 (8 pp.), http://dx.doi.org/10.1155/2008/ 358748 14. Xia G, Tawhai M, Hoffman E, et al. Airway wall stiffening increases peak wall shear stress: a fluid–structure interaction study in rigid and compliant airways. Ann Biomed Eng 2010; 38(5): 1836–1853. 15. Philips CG and Kaye SR. On the asymmetry of bifurcations in the bronchial tree. Respir Physiol 1997; 107: 85–98. 16. Nakano Y, Wong JC, de Jong PA, et al. The prediction of small airway dimensions using computed tomography. Am J Respir Crit Care Med 2005; 171: 142–146. 17. Durst F, Ray S, Unsal B, et al. The development lengths of laminar pipe and channel flows. J Fluid Eng: T ASME 2005; 127: 1154–1160. 18. Horsfield K, Dart G, Olsen DE, et al. Models of the human bronchial tree. J Appl Physiol 1971; 31: 207–217. 19. Calay RK, Kurujareon J and Holdø AE. Numerical simulation of respiratory flow patterns within a human lung. Resp Physiol Neurobiol 2002; 130: 201–221. 20. Tgavalekos NT, Tawhai M, Scott Harris R, et al. Identifying airways responsible for heterogeneous ventilation and mechanical dysfunction in asthma: an image functional modeling approach. J Appl Physiol 2005; 99: 2388–2397. 21. Mauroy B, Filoche M, Andrade J, et al. Interplay between geometry and flow distribution in an airway tree. Phys Rev Lett 2003; 90(14): 148101. 22. Tarran R, Button B, Picher M, et al. Normal and cystic fibrosis airway surface liquid homeostasis – the effects of phasic shear stress and viral infections. J Biol Chem 2005; 280(42): 35751–35759. 23. Sidhaye VK, Schweitzer KS, Caterina MJ, et al. Shear stress regulates aquaporin-5 and airway epithelial barrier function. Proc Natl Acad Sci USA 2008; 105(9): 3345–3350. 24. Evans DJ and Green M. Small airways, a time to revisit? Thorax 1998; 53: 629–630. 25. Button B, Picher M and Boucher RC. Differential effects of cyclic and constant stress on ATP release and
mucociliary transport by human airway epithelia. J Physiol 2007; 580(2): 577–592.
Appendix 1 The CFD model A CFD model has been chosen to examine velocity profiles and shear stress distributions in a model of small airways. The CFD technique approximates and numerically solves the incompressible Navier–Stokes fluid flow equations (1) and (2). There is no analytical solution for these equations in a flow system due to the complexity of that geometry (such as a bifurcating network of tubes). Therefore, the equations have to be solved numerically using a computer software program, and such techniques have been developed and refined over the last 40 years r u=0 ru ru = rp + mr2 u
ð1Þ ð2Þ
where u is the velocity vector (m/s), p is the pressure (Pa), r is the air density (kg/m3) and m is the air molecular viscosity (Pa s).
Wall shear stress The wall shear stress t only occurs where flow interacts with a surface (i.e. an airway wall), and it is a derived quantity defined as the product of molecular viscosity with the local perpendicular velocity gradient, t = m(du/dy). Poiseuille flow assumes ideal conditions of a long straight pipe18 of constant circular diameter and a steady, laminar flow where the velocity profile has developed to a parabolic velocity shape; velocity is inversely proportional to the square of the pipe diameter. The wall shear stress may be calculated5,14 as t=
32 mQ pd3
where Q is the airway volumetric flow rate (m3/s) and d is the airway diameter (m), and thus, the wall shear stress is inversely proportional to the cube of the pipe diameter. The degree to which conditions depart from the ideal, as in a flow bifurcation with pipe lengths of very few diameters, may be readily compared by nondimensionalising wall shear stresses to those derived from Poiseuille.
Downloaded from pih.sagepub.com at University of Bath - The Library on June 18, 2015
