Flow Characterization in Healthy Airways and Airways with Chronic Obstructive Pulmonary Disease (COPD) during Different Inhalation Conditions

The flow of air in the central and lower generations of the human airways results in complex secondary flows, which heavily affect the deposition of aerosol such as PM2.5 in the human airways. To figure out the flow dynamics in healthy airways and airways with chronic obstructive pulmonary disease (COPD), the flow distributions were predicted numerically. In a bronchial tube geometry (G5-G8), COPD was represented by an axisymmetric constriction in G6-2 of a fourgeneration airway. Real and unsteady breathing inhalation curves under rest, light activity, and moderate exercise were applied in the simulations. The velocity distribution, recirculation, and stagnation zones in the generation G6-2 affected by COPD together with the effects of sudden pressure drop in the constriction were investigated. Two symmetrical vortices (Dean vortices) were formed at the cross-sections of the geometry affected by COPD and their intensities were dependent on the Reynolds number. The mas flow rates in the healthy generation (G6-1) went to as high as 3.10 times the mass flow rate in the generation affected by COPD (G6-2). The simulations showed that inflammation of the bronchus greatly affected pressure distribution and the rates of mass flow in the human airway during inhalation. The overall pressure drop in the bifurcation affected by COPD was 3.38 Pa for rest, 10.91 Pa for light activity and 21.39 Pa for moderate exercise. A combination of these factors make COPD patients very susceptible to adverse health effects caused by inhalation of PM2.5.


A
Amplitude of the real breathing curve (Reynolds number) D Diameter (mm) F Force (N) G Generation of Weibel's airway P Pressure (Pa) Q Flow rate (L s -1 ) R Radius (mm) U i Airflow velocity (m s -1 ) x i Characterizes the Cartesian coordinates σ ij Shear stress (Pa) Re Reynolds number t Time (s) U mean Average velocity at the airway's entrance (m s -1 )

INTRODUCTION
Global initiative for Chronic Lung disease (GOLD) defines chronic obstructive pulmonary disease (COPD) as a human airway condition associated with persistent and limited airflow which is usually progressive and associated with higher inflammatory response to toxic particles and gases (Rabe et al., 2007). This condition is currently widespread, as indicated by the data from global burden of diseases in the year 2015, where 5% of all mortality (3.17 million) in the world was due to COPD (Mathers and Loncar, 2006). Epidemiological reports from World Health Organization (WHO) further estimated that there were 251 million cases of COPD in the year 2016 (Mathers and Loncar, 2006). Currently, COPD is ranked as the fourth leading cause of death in the world. Projections showed that it will be the third leading cause of death by 2030 (Adeloye et al., 2015). Projections further showed that the number of people suffering from COPD are likely to reach 7.8 million in the year 2030 (Mathers and Loncar, 2006). This condition has been linked to tobacco smoking as well as fine particulate (PM 2.5 ) and other toxic pollutants from both bio-genic sources and non-biogenic sources for instance combustion processes (Pipal et al., 2014;Zhu et al., 2017). COPD leads to reduced quality of life, burden in terms of cost associated with the management of the disease and also death (Li et al., 2018).
As noted by Kolanjiyil and Kleinstreuer (2017), the air dynamics inside the human airways will determine the deposition patterns of pollutants, especially fine particulates PM 2.5 which are usually suspended in the ambient air Zhu et al., 2017;Wang et al., 2018). The emergence of biotechnology industrial activities, e.g. waste recycling, detergent as well as food industries in the recent years has increased human exposure to bio-aerosols. Information of the deposition patterns and efficiencies for the bio-aerosols can provide insight in carrying out epidemiological investigations (Douwes et al., 2003). Several studies have been made to estimate the deposition efficiencies of particulate matter in the human airways, as presented in Table 1 (Luo et al., 2007;Feng and Kleinstreuer, 2014;Chen et al., 2018). In vivo lung flow measurements in both healthy and defective human airways are hard to administer due to the complexity of both the human airways and the breathing process. The application of in vitro flow measurements for both healthy individuals and those affected by COPD is limited to the upper airways (mouth to G3). The limitation is due to the clarity of scanned images of the human airways and reproducibility of the scanned image which are very small in dimensions after G3. As such, there is a need to be innovative in order to estimate the flow rates in the peripheral human airways (Delvadia et al., 2012).
Computational fluid dynamics (CFD) has been widely applied by researchers to predict physical phenomenon in the natural environment (Chen, 2001a, b;Chen et al., 2011;Panagopoulos et al., 2011). Therefore, the limitations of flow measurements in the human lungs discussed earlier can be solved through the application of CFD (Chen et al., 2014;Kolanjiyil and Kleinstreuer, 2017). A study conducted by Weibel (1963a) provided geometric data for the human airways as symmetrical bifurcations with 23 generations. A number of studies have been published on the application of this geometry in flow models for the simulation of flow effects in healthy airways and those which have been deformed by asthma, COPD, or tumors (Kim and IGLESIAS, 1989;Zhang and Papadakis, 2010;Chen et al., 2018).
Lack of symmetry in buckled human airways has been shown to cause massive changes in terms of flow rates and secondary flow patterns (Liu et al., 2015). In this study, flow velocity distribution in G5-G8 of the Weibel's bifurcation for a normal airway and a COPD case were studied. Yang (2006), applied steady flow distribution for a simulation during an earlier study on flow velocity distribution in human airways affected by COPD. The results showed that maximum velocity, recirculation, and stagnation zones were higher in an obstructed airway compared to a normal airway because of the secondary flow velocities caused by the constriction (Yang et al., 2006). Previously, it has been shown that airflow pattern had an impact on the ease of breathing and deposition patterns for particulate matters in the human lungs (Huang and Zhang, 2011;Feng and Kleinstreuer, 2014). Studies have been conducted about the velocity distribution during a light activity using an unsteady velocity distribution at the inlet (Chen et al., 2012).
COPD patients are faced with reduced life quality and shortened life expectancy, mainly due to obstructions which alter the airflow in the human airways, hence making it difficult to breathe normally (Luo et al., 2007). Consequently, prediction of airflow in the human airways affected by COPD is vital in order to assess the effects of the inflammation on the mass flow rates of air inside the airways (Yang et al., 2006). Chen et al. (2012) studied particle deposition in a COPD case during a light activity. However, studies on airflow patterns for human airway bifurcation affected by COPD under real breathing conditions for rest and moderate exercise conditions have not been considered yet. Therefore, as an improvement of the study by Chen et al. (2012), airflow simulations were performed for COPD and healthy airways under real transient flows during rest, light activity, and moderate exercise in the present study. The three breathing conditions were investigated and the properties of diverging patterns of airflow and their velocities of the airways were evaluated using numerical methods. The results from this study enable us to understand the effects of different activity levels on the mass flow rates through the distorted airways such as for individuals affected by COPD.

Geometry of the Model
Currently, there are two major approaches applied in the construction of lung models. The first approach involves digitally reconstructing the human airway geometry from detailed imaging data (Katz et al., 1999). Working on computed tomography (CT) scanned data will definitely offer more realistic results (Comer et al., 2000). The major drawback associated with this method is the limited prospect of reproducing scan-based models because the imaging procedures are both costly and time consuming (Kolanjiyil and Kleinstreuer, 2017). Additionally, obtaining scan based models is presently limited to the mouth cavity and a few upper airway generations (G0-G3) due to limited resolution of the scanned images for central and lower generation of the human airways (Kolanjiyil and Kleinstreuer, 2017).
The second approach involves basing the human airway geometry on structures that have satisfied certain modelling requirements (Chen et al., 2018). The most commonly used geometries for the modelling airflow and particulate matter deposition in the human airway are Weibel's model and Horsfield's model (Weibel, 1963b;Horsfield et al., 1971). For this simulation, the realistic bifurcation based on the Weibel's model was preferred due to its accurate numerical solutions and efficiency in terms of computational time and resources (Piglione et al., 2012). Two sets of human airway models were developed in order to validate the approach used for this study and investigate the airflow characteristics (Chen et al., 2012;Sul et al., 2014). A healthy geometric model of G5-G8 of the human airways, and a COPD geometry (G5-G8) wherein G6-2 was affected by COPD were developed. The generations were selected because there was a need to compare the velocity distributions with the results from Chen et al. (2012) for validation purposes, as well as avoid the turbulent flow experienced in the larynx and upper section of the airway (Katz et al., 1999). The cylindrical models of G5-G8 with a bifurcation angle of 70° is shown in Fig. 1 (Weibel, 1963a). The meshes for the geometry affected by COPD and the healthy geometry are as shown in Figs. 2(b) and 2(c) respectively.
The geometry of a human airway affected by COPD can be represented by a tube with a constriction along its length to represent the region which has been affected by the inflammation (Yang et al., 2006;Chen et al., 2012). The constriction reduces the cross-sectional area of the airways, thereby increasing the resistance of airflow in the affected airway path. The geometry applied for this simulation was similar to that used by Chen et al. (2012) to simulate deposition of 5 μm diameter particles in the human airways. Finite element analysis was applied to investigate and compare the airflow patterns in a healthy airway geometry and a COPD case (Feng and Kleinstreuer, 2014;Chen et al., 2018). The parameters of the geometries used for the simulation are summarized in Table 2. The geometry of the generation affected by COPD is defined such that constriction and subsequent expansion are symmetrical about the axis of the airway. The axial diameter of constriction is half of its original size. The openings at the end of G6-2 are smoothly connected with the constriction using an arc. The other respective parameters of the model geometry are provided in Table 2.

Governing Equations and Boundary Conditions
This study focused on the gas flow fields only during inhalation in healthy airways and airways with COPD. It was assumed the flows in the airways were isothermal and laminar, this implies, in turn, the flows are incompressible. The airways were assumed to be smooth rigid tubes with no-slip on the tube surfaces (Yang et al., 2006;Chen et al., 2012;Li et al., 2013). The density (ρ) and kinematic viscosity (μ) of the air inhaled were all assumed to be constant. The basic conservation equations of continuity and momentum (Navier-Stokes) equations for the flow of a gas through the geometry are expressed in Fig. 1 (Inthavong et al., 2010;Rahimi-Gorji et al., 2015).
Different statuses of breathing depend on an individual's level of activity; they include rest, light activity, and moderate exercise . Real breathing curves from the studies of  and Chen et al. (2012) on transient inhalation were used in this study. Therefore, the variation of mean velocity of air at the    Fig. 3 Chen et al., 2012). The three intensities of activities and the respective factors for inhalation settings are shown in Table 2. The mean velocity at the inlet of the parent tube of the bifurcation was also described as shown below (Li et al., 2018) where D represents the diameter of the bronchus at the inlet, G stands for the generation of the airway according to Weibel's model, and Q in is the airflow rate at the entrance. For the numerical analysis, a completely-formed parabolic velocity distribution was applied at the entrance of the geometry. ANSYS 16.0 provides the user-defined function (UDF) which was used to create the velocity distribution profile at the inlet of the geometry. The real velocity distributions at the inlet of the parent tube with time for the three inhalation statuses were curve-fitted based on the data in the real inhalation curves during the respective inhalation statuses at G5 Chen et al., 2012) and presented in Table 3. The coefficients of determination of the curve fittings are greater than 0.997 (i.e., R 2 > 0.997), showing excellent fitting of the curves. This was an improvement of previous studies where a sinewave was used to represent the variation of velocity at the inlet of the airways with time (Chen et al., 2018).

Reynolds Number
To define the velocity at the inlet of the bifurcation, the Reynolds number of air at the entrance of the geometry (Re d ) was specified and applied for the numerical simulation . The Reynolds number is given as: where ρ is air density (= 1.225 kg m -3 ), U mean is the average velocity at the entrance of G5, μ is the dynamic viscosity (= 1.789 × 10 -5 kg m -1 s -1 ), and D represents the size of the diameter at the entrance of G5 (=3.5 × 10 -3 m) (Weibel, 1963b). At the mouth region, which was the entrance to the human airways, air velocity was uniformly distributed. However, when the flow proceeded from G0 into other generations, the velocity profile became a parabolic distribution due to drag force developed on the walls of the airways where the no-slip boundary condition was imposed (Zhang and Papadakis, 2010). The airflow velocity profile at the entrance of G5 had a parabolic distribution with the maximum flow velocity (U max ) being at the middle of the cross-sectional area of the airway and was double the mean flow velocity of the airway at the entrance of the airway (U mean ). For a fluid whose density and dynamic viscosity were fixed, the average and maximum flow velocities (U mean end U max ) varied proportionally with the Reynolds number (Re) at the entrance of G5. The assumption was made due to its practicability and also for comparability with other studies Zhang and Papadakis, 2010;Chen et al., 2012)

Numerical Method and Simulation
A three-dimensional four-bifurcation geometry was created using the graphical software SOLID WORKS. The mesh for the geometry was generated based on the tetrahedral elements (unstructured grid system) using the commercial software ANSYS 16.0. The physical phenomena were predicted by solving the equations of fluid flow. The equations were discretized using an inbuilt finite-volume method, and the simulations of the single-phase flow in the geometry with four generations followed the SIMPLEC algorithm. SIMPLEC algorithm has gone through a significant improvement over time and is currently capable of handling a low Reynolds number flow with an acceptable accuracy. A three-order MUSCL scheme was used to discretize the momentum of the air, whereas the gradient selection was established according to Green-Gauss Node option. The criteria for convergence was based on pressurevelocity coupling. In this study, the relaxation coefficients for momentum and pressure were set to be 0.75.
During the simulation, rest had 60 time steps, light activity had 57 time steps, while moderate activity had 34 time steps. To ensure the smoothness of the positive half real inhalation cycle, the time step of 0.035 s (Δt) was carefully chosen. The solutions for the flow fields at each time step were iterated or calculated repeatedly until they converged. At present, it was presumed that the calculation was convergent when the mass residue left in the continuity equation was less than 10 -5 (Chen et al., 2018).

Model Validation and Grid Independence
Numerical validation is an important step in the simulation. Therefore, in order to validate the developed numerical model, the results from the numerical simulations in this study were compared with those of Chen., (2012). Two grid systems with tetrahedral meshes were built for the healthy geometry and the geometry affected by COPD, as shown in Figs. 2(b) and 2(c). The computational area included four generations (G5-G8), and the average Reynolds number at the entrance of G5 was 164.3 for rest, 362.4 for light activity, and 606.4 for moderate exercise, as shown in Table 4 Chen et al., 2012).
The velocity contours during a light activity for G5 to G7 at an inhalation time of 1.0 s during a light activity are shown in Fig. 4(a), while Fig. 4(b) shows the pressure contours at the generations affected by COPD. Overall, the velocity contours and pressure contours in the present study were in line with the results of Chen et al. (2012), thereby validating the present simulations. To ensure the accurate simulations, the grid independence was also tested by comparing the mass flow rate and mean pressure at three different cross-sections in G6-2 (airway with COPD) and  G6-1 (healthy airway) (Fig. 5(a)) under three different grid systems; they are Mesh 1 with 1,028,518 elements, Mesh 2 with 1,500,518 elements, and Mesh 3 with 2,043,737 elements. The study was based on a real breathing at the inhalation time of 1.0 s. The predicted mass flow rates ( Fig. 5(b)) and mean pressures (Fig. 5(c)) depict that the values of mass flow rate and mean pressure drop based on the three grid systems are very close to each other. This implies, in turn, that Mesh 2 with 1,500,518 tetrahedral grids was enough for the simulations based on the geometry with COPD and thus employed for the simulations of flow phenomena in the airways. While the healthy geometry had 1,312057 elements and 534,266 nodes, the geometry with COPD had 1,500,186 elements and 278,270 nodes in order to capture the complex fluid flow phenomena at the generation affected by COPD. The two configurations were carefully chosen for both healthy and COPD case in order to guarantee accurate and reliable results.

Mass Flow Ratios and Maximum Velocities in Airways
In order to figure out the flow behaviors in healthy airways and airways with COPD, air mass flow rates and the maximum velocities in G5-G7 under three different breathing situations are given in Table 5. For the healthy airways, the flow rate in G6-1 and G6-2 are equivalent, due to the symmetrical geometry. The maximum velocity for the healthy airways always develops at G5 because the total cross-sectional area was the lowest at that generation. It is noteworthy that the mass flow rate and the maximum velocity in G7-2 or G7-3 are greater than those in G7-1 or G7-4, as a consequence of centrifugal force (Luo et al., 2007). For the airways with COPD, G6-2 with geometrical contraction has a lower mass flow rate when compared to G6-1, irrespective of the breathing status. Unlike in healthy airways, the maximum velocity in airways with COPD This study Chen et al., 2012 (a) (b)

Fig. 5.
A schematic of (a) block structure geometry for G5-G7 wherein G6-2 is affected by COPD, and (b) mean pressure, and (c) mean mass flow profiles at different cross-sections under three different grid systems at light activity and 1.0 s. appears in G6 rather than G5. The mechanism causing different location of the maximum velocity will be discussed later.
The profiles of mass flow ratio in healthy G7 and G7 with COPD are shown in Fig. 6. The mass flow ratio is defined as: where ṁ i is the i th branch mass flow rate in G7, and ṁ e is the total mass flow rate at the entrance of G7. For the healthy airways under rest, the difference in the mass flow ratio is not significant. However, the higher the respiration rate for inhalation (Table 4), the more pronounced the difference. The mass flow ratios of G7-2 and G7-3 increase with increase in the Reynolds number, whereas those of G7-1 and G7-4 decrease. The ratios in G7-2 and G7-3 are of 0.258 during rest, 0.269 during light activity, and 0.272 during moderate exercise. For the COPD model, the highest mass flow ratio of 0.398 is exhibited in G7-2 during light activity, while the lowest mass flow ratio of 0.073 is observed in G7-4 for moderate exercise. For G7-1 and G7-3, there is an increasing trend in the ratio with increasing Reynolds number, while for G7-4 the ratio is indirectly proportional to the Reynolds number as shown in Fig. 6. These differences in the mass  flow ratios in the COPD are due to the different levels of resistance of air-flow induced by the presence of the constriction, as discussed in a later section.

Airflow in the Healthy Airways
To figure out the influence of COPD on the flow dynamics in the airways, the mid-plane velocity distributions in the healthy airways at three different inhalation conditions, namely, rest, light activity, and moderate exercise, are shown in Fig. 7. For the three inhalation statuses, the maximum airflow velocities in the sinusoidal waveform in Fig. 3 were used to show the flow velocity distributions . It can be seen that corresponding to the statuses of rest, light activity, and moderate exercise, the inhalation times are approximately 1.0, 1.0, and 0.6 s, respectively.
For the healthy airways, the presence of curvatures due to the bifurcations led to the concentration of streamlines of the flow towards the inner walls, and thus higher mass flow rates at these zones are observed shown in Fig. 7. This effect is due to the induced centrifugal forces in the geometry of airways (Zhang et al., 2001Kim and Lee, 2009). The higher the Reynolds number, the more significant the curvature effect or the centrifugal force on the flow concentration. This is the reason why the low velocity zones close to the outer walls in G6 expanded markedly when the mass flow rate increased. Based on the centerline of G5, the geometry of the healthy airways was symmetrical, so it was not surprising that the airflows were also symmetrical, regardless of the inhalation status. However, the flow field in G6-1 or G6-2 was not symmetrical, as the consequence of the centrifugal forces. Hence, the airflows in G7 were not symmetrical either, even though the geometry along the centerline of G6-1 was symmetrical. This result was in good agreement with the findings of Zhang and Papadakis (2010). The curvatures in the flow paths after the first bifurcation (in this case between G5 and G6), the fragmentation of flow at the bifurcations, and the momentum of the airflow led to skewed velocity distributions. The lowest velocities were at the outer edge of the bifurcation in the daughter branches, even when the velocity distribution in the mother branch was symmetrical (Guha et al., 2016). The lowest mass flow rates in each generation were located in the outermost airways (i.e., G7-1, G7-4, G8-1 and G8-8 in Fig. 2(a)), whereas the highest mass flow rates were exhibited in the innermost airways (i.e., G7-2, G7-3, G8-4 and G8-5 in Fig. 2(a)) for all the inhalation statuses. This phenomena was markedly amplified with an increase in Reynolds number, resulting in even more skewed mass flow ratios at a higher Reynold number.

Airflow in Airways with COPD
For the airway with COPD, its geometry is a contraction for flow in nature. This implies, in turn, that the air in the channel has a higher velocity accompanied by a higher head loss due to friction and geometry change. The midplane velocity distributions in the airways with COPD at the three different inhalation statuses are displayed in Fig. 8. Overall, it can be seen the presence of stagnation and recirculation zones in G6-2, namely, Z1, Z2, and Z3, due to the constriction, and an increase in Reynolds number enlarged the sizes of the stagnation and recirculation zones. So the zones were largest in the case of moderate exercise (Fig. 8(c)) (Yang et al., 2006;Luo et al., 2007). Fig. 8 depicts that the location of the maximum velocity was sensitive to the inhalation status. For rest and light activity, the maximum velocities were in the airway with COPD (G6-2), and their values were 2.49 ( Fig. 8(a)) and 4.62 m s -1 (Fig. 8(b)), respectively. On the other hand, the maximum velocity in the case of moderate exercise was 7.86 m s -1 , which appeared in the healthy airway G6-1 (Fig. 8(c)). These phenomena could be explained in terms of energy loss and flow change in the system, as a result of the gradual contraction and expansion of the airway affected by COPD.
The energy loss due to the geometry change in G6-2 was by far smaller when compared to the loss due to friction in an entire generation; hence, the former pertained to the minor loss, whereas the latter belonged to the major loss. In G6-2, the minor loss (h COPD ) was caused by a gradual contraction followed by a gradual expansion, and was expressed by where K is the resistance coefficient, v is the inlet average velocity, g is the gravitational constant, and the subscripts gc and ge stand for gradual contraction and gradual expansion, respectively. For inhalation during rest, the minor loss was small, stemming from the low velocity of airflow. Under such a situation, the maximum velocity observed in the cross-section B'-B was mainly due to the contraction of the cross-sectional area which increased the local flow ( Fig. 8(a)). In contrast, for moderate exercise, the inlet velocity was large so the minor loss became important. This gave rise to more airflow into the healthy airway. This was responsible for the maximum velocity developed in G6-1 (Fig. 8(c)). In regard to light activity, the minor loss and geometry change almost played an equivalent role in determining the maximum velocity. Consequently, the maximum velocities in G6-1 and G6-2 were close to each other. Nevertheless, the maximum velocity appearing in G6-2 revealed that the factor of minor loss was slightly higher than that of geometry change.

Velocity Vectors and Contours in the Airway with COPD
In order to understand the effect of COPD on the airflow in the airways, the velocity vectors and the axial velocities in G6-2 at three different cross-sections A'-A, B'-B, and C'-C under the three inhalation statuses are shown in Fig. 9. The cross-sectional areas of A'-A, B'-B, and C'-C were 6.16, 3.07, and 6.16 mm 2 , respectively. In A'-A, the regions of high velocity were closer to side A (i.e., the inner wall in G6-2), irrespective of the inhalation status. This is attributed to the centrifugal force induced after the bifurcation (Yang et al., 2006;Luo et al., 2007). Meanwhile, it is noted that the secondary velocity vectors resulted into two symmetrical vortices, termed the Dean vortices, developed in the crosssections due to skewed flow velocities in G6-2 (Dean, 1927). The flow rate during rest was lowest ( Fig. 9(a)), so the vortices were less pronounced compared to those formed during the other two breathing conditions (Figs. 9(b) and 9(c)). Accordingly, it is recognized that the development of Dean vortices is strongly dependent on flow velocity or Reynolds number, which agreed well with the findings of Kim and Lee (2009).
In G6-2, the centrifugal force was highest at A'-A. As a result, the dean vortices were present for all inhalation statues (Fig. 9). The highest flow velocities at the constriction were observed during moderate exercise ( Fig. 9(c)) and this finding was in line with the findings by Choi et al. (2007) and Chen et al. (2012). On account of high resistance of flow in the constriction (i.e., B'-B), the dean vortices were not seen during rest, but was evident under light activity and moderate exercise. Unlike in A'-A, the axial velocity contours in B'-B were located at the center of the crosssection. This revealed that the impact of the centrifugal force on the flow structure in B'-B was not as significant as in A'-A. The velocity vectors of the secondary flow were different before and after the significant reduction in cross-sectional area, resulting from the variation in flow velocity and centrifugal forces (Guha et al., 2016). Convergence of the airflow due to the contraction generated higher flow velocities for all inhalation conditions. Centrifugal forces were weakest in C'-C. This was ascribed to two factors: (1) the greater distance from the bifurcation at the end of G5 and (2) the energy loss by the contraction. Consequently, the dean vortices only present for moderate exercise (Fig. 9(c)). This was consistent with the study of Sul et al., (2014) in which the airflow resistance was used to characterize obstructive lung diseases. Meanwhile, the subsequent increase in the tube's radius led to reduced flow velocity in C'-C. This change in the geometric shape (a) resulted into complex secondary flow patterns in the crosssection (Zhang et al., 2001). The flow velocities were highest during moderate exercise, which gave rise to both the velocity vectors and axial velocity contours showing the vertically symmetrical vortices.

Variations in Average Velocities and Pressure Drop in the Generation Affected by COPD
Due to the constriction, there was an accelerated flow velocity phenomenon which led to high average velocities at the center of the constriction (cross-section B'-B). Subsequent expansion of the geometry to its original crosssectional area led to a decrease in the average velocities to their original values as shown in Fig. 10(a). The average velocity of 4.87 m s -1 was highest at the center of the constriction (cross-section B'-B) during moderate exercise. Pressure drops in the human airways are attributed to frictional and minor losses. Frictional losses, which account for the biggest share in pressure drop, usually result from drag as air flows along the length of the generation. Alternatively, minor losses occur at the bifurcations and at the constriction (for a COPD case) due to the gradual change in the geometry of the airways (Kang et al., 2011). Based on the principle of conservation of energy, the pressure drop across a single generation is composed of three items: (1) instantaneous acceleration of the airflow within the control volume; (2) energy loss due to viscosity of air; and (3) the head losses in the airway due to friction , and is expressed as where Ω is control volume of the geometry under study, ∆P is the average static pressure drop, ∆P v is the viscous pressure drop, and ∆P k is the change in pressure caused by change in kinetic energy. The gradual reduction in crosssectional area led to low mean pressure and high average velocity in the constriction, as indicated in Fig. 10. Pressure remained low after the subsequent expansion of the airway. This is similar to an observation made by Qi et al., (2014), whereby the low pressure regions coupled with high shear stress were observed in bifurcations affected by bronchial stenosis. The profiles of the mean pressure for a COPD case in the three cross-sections A'-A, B'-B and C'-C are shown In Fig. 10(b). Overall, the higher the Reynolds number, the larger the pressure drop observed. During rest (Re = 164.3), the pressure drop between A'-A and B'-B was 3.38 Pa for a COPD case as compared to 0.43 Pa for the healthy geometry. For light activity (Re = 362.4), the pressure drop was 10.91 Pa for a COPD geometry as compared to 1.10 Pa for the healthy geometry. Lastly, during moderate exercise (Re = 606.4), the pressure drop was 21.39 Pa for a COPD case as compared to 2.60 Pa for the healthy geometry. This characteristic was also demonstrated by Jain et al. (2013) where the study proved that the presence of tumors in the human airways would lead to a sudden pressure drop, thereby reducing the removal of bacteria from the airways resulting in weakened defenses. Patients affected by COPD suffer from impaired bronchial clearance due to low mass flow rates and high pressure drops caused by the inflammation. As a consequence, their airways suffer from increased deposition and accumulation of particulate matters among other toxins found suspended in the air . For this reason, particulate matters suspended in the atmosphere affect individuals living with the COPD condition in a more adverse way compared to the healthy population (Yang et al., 2006;Luo et al., 2007). The present (a) (b) Fig. 10. Profiles of (a) average velocity and (b) mean pressure in cross-sections A'-A, B'-B, and C'-C for the COPD case at the peak flows for different inhalation conditions. study has provided useful insights into the recognition of the fluid dynamics in healthy airways and airways with COPD under different inhalation conditions.

CONCLUSIONS
Computational fluid dynamics has been successfully applied to evaluate the airflow fields in healthy airways and airways with COPD where four generations, G5-G8 under three different inhalation statuses, namely; rest, light activity, and moderate exercise, were considered. Dean vortices were observed in the area next to the inlet of G6 and were nearer to the inner walls, for both the healthy and the COPD case. For the healthy airways, the mass flow rates were equivalent for corresponding generations in the geometry as compared to the geometry with COPD where higher mass flow rates were witnessed in the healthy generation compared to the generation affected by COPD. Due to low minor losses during rest, the maximum velocity was exhibited at the constriction. The balancing of the effects of minor losses and geometry change during a light activity led to the maximum velocities in the generation with COPD and the healthy generation being almost equal. The minor losses became significant for moderate exercise, hence the highest velocity was located in the healthy airway. The presence of the constriction in G6-2 strongly alters the airflow patterns due to variations in minor losses and hence the resistance of flow. As a result, the mass flow rates in G6-1 were 2.78, 3.04, and 3.10 times that in G6-2 for rest, light activity, and moderate exercise, respectively. For the airways with COPD, pressure drop was highest at the location of the constriction. The highest pressure drop was during moderate exercise in the geometry affected by COPD at 21.39 Pa as compared to 2.60 Pa for a healthy airway. The intensity of Dean vortices and the degree of asymmetry of flow in the geometry affected by COPD, were more substantial during light activity and moderate exercise as compared to rest. These findings have important implications on the effects of activity level and presence of a constriction on air transportation within the airways affected by COPD.