Abstract

Better understanding of stresses and flow characteristics in the human airways is very important for many clinical applications such as aerosol drug therapy, inhalation toxicology, and airway remodeling process. The bifurcation geometry of airway generations 3 to 5 based on the ICRP tracheobronchial model was chosen to analyze the flow characteristics and stresses during inhalation. A computational model was developed to investigate the airway tissue flexibility effect on stresses and flow characteristics in the airways. The finite-element method with the fluid-structure interaction analysis was employed to investigate the transient responses of the flow characteristics and stresses in the airways during inhalation. The simulation results showed that tissue flexibility affected the maximum airflow velocity, airway pressure, and wall shear stress about 2%, 7%, and 6%, respectively. The simulation results also showed that the differences between the orthotropic and isotropic material models on the airway stresses were in the ranges of 25–52%. The results from the present study suggest that it is very important to incorporate the orthotropic tissue properties into a computational model for studying flow characteristics and stresses in the airways.

1. Introduction

Flow characteristics in the human respiratory airways are very important for studying particle transportation in many applications ranging from aerosol drug therapy to inhalation toxicology as well as gas exchange in the airways. A better understanding of the stresses and strains and theirs influences on our airways is also important, since the airwaysare shaped by the complex mechanical environment even in the uterus. This complex mechanical environment continues to influence and alter the mature airways in healthy and diseased people [1]. Many in vivo and in vitro models have been developed to study the effect of the mechanical stresses on the airways. The results from these models showed that the high peak airway pressure caused airway distention and an increase in the level of an inflammatory mediator, cytokine [2, 3]. In addition, high airway pressure can contribute to an increase in the thickness of the airway epithelial layer [4, 5]. Since it is very difficult to experimentally measure the flow characteristics (airflow velocity, airway pressure, and wall shear stress (WSS)) and airway stresses in the real environment, the computational model can provide useful information about flow characteristics and stresses in the airways.

Many researchers have developed computational models to investigate airflow and particle deposition in the airways [612], including the effect of airway diseases, such as tumors [610], asthma [11], stenosis [13], and chronic obstructive pulmonary disease (COPD) [12, 14]. The effects of carinal shape [15] and cartilage rings [16] were also studied. However, these models were developed assuming that the airways were rigid and could not be deformed or distended by fluid force from the airflow. In general, the airways are flexible and are composed of many soft tissues. They can also be distended by airflow during inhalation [17]. The literature review indicates that there is no computational model that considers the effect of airway tissue flexibility on the flow and tissue characteristics.

In the present study, a fluid-structure interaction (FSI) analysis along with the finite-element method was employed to investigate the flow characteristics and stress distributions in the airways during inhalation. Both orthotropic and isotropic material models were used to represent airway flexibility to study stress distributions in the airways. The results between the fluid analysis and FSI analysis are also compared.

2. Materials and Methods

The transient interactions between airflow and airways during inhalation were investigated solving two coupled sets of the governing equations with associated boundary conditions. The governing equations for airflow and airways are briefly described below.

2.1. Governing Equations for Airflow

The governing equations for transient airflow are Navier-Stokes equations on a moving mesh with the assumption of incompressible flow. These equations govern the principles of mass and momentum conservation and are described below using Einstein’s repeated index convention [18].

Conservation of mass Conservation of momentum In these equations, represents the moving mesh location, is the metric tensor determinate of the transformation, that is, the local computational control-volume size, is fluid density, is fluid pressure, is fluid viscosity, and is fluid velocity.

2.2. Governing Equations for Airways

The governing equations for movement of the airways during inhalation are the time-dependent structural equations and are described below using Einstein’s repeated index convention [19].

Equation of motion Constitutive relations In the equation above, is the stress in each direction, is the body force, is density, and is the displacement, is the elasticity tensor, and is the strain in each direction.

2.3. Computational Method

The effect of fluid pressure on a structure is significant, especially if the structure is flexible such as human airways. The finite-element software ANSYS [20] was employed to solve this fluid-structure interaction problem. The solution to the fluid-structure interaction problem can be obtained by solving the governing equations for airflow and airways consecutively. At each time step, the algorithm begins by solving the airflow equations to obtain fluid pressure. Airway equations are then solved for the displacement using the fluid pressure as an external force. The airflow equations are then solved again to obtain the fluid pressure after the airway displacement changes the fluid boundaries. This loop continues until both fluid pressure and airway displacement converge for each time period (see Figure 1).

2.4. Computational Domains

This study focused mainly on airway generations 3 to 5 for two reasons. First, these airway generations have less cartilage plates and no rings when they are compared to the proximal generations; therefore the airways can be assumed to be smooth [21]. Second, the diameters of these airways do not change as a function of a lung volume but depend on a transmural pressure across the airways [17]. The geometric dimensions of airway generations 3 to 5 (see Figure 2) used in this study were based on the ICRP [22] tracheobronchial geometry, and airway thickness for each generation was based on measurements by Habib et al. [23]. The branching angle of the bifurcation was 70° based on the morphological data of Horsfield and Cumming [24]. The corresponding geometric diameter, length, and thickness of the bifurcation are given in Table 1. The surface geometry of the model was constructed based on the physiologically realistic bifurcation (PRB) model suggested by Heistracher and Hofmann [25]. This double bifurcation geometry was previously implemented in a study by Longest and Vinchurkar [26], which evaluated the effect of transitional and turbulent flow on particle deposition in rigid airways.

2.5. Computational Models and Boundary Conditions

The computational domains of the bifurcation were created in ANSYS [20]. Due to symmetry, only one half of the domains was constructed. The solid domain was the airways with finite thickness, and the fluid domain was the internal volume of air in the airways. Solid elements, BRICK45 [20], were used to represent the solid domain, and fluid elements, FLUID142 [20], were used to represent the fluid domain. A structural hexahedral mesh was employed to provide a high-quality airflow solution, as suggested by Longest and Vinchurkar [27] and Vinchurkar and Longest [28]. A mesh-independence study was performed on both solid and fluid domains to confirm that a fine enough element had been used to represent the solid and fluid domains. Changes in maximum pressure and velocity were used as convergence criteria for the fluid domain, and changes in maximum displacement and von Mises stress were used as convergence criteria for the solid domain. A converged model was obtained when changes in those criteria were less than 4%. Having performed the mesh-independence study, the airflow velocity from the finite-element model in the airway generation 4 was then compared with the experiment by Zhao and Lieber [29]. Good agreement was obtained between the simulation and experimental results [30]. Figure 3 shows the converged mesh that was further used for the analysis of predicting stress distributions.

The inlet boundary condition of the fluid domain was a waveform during normal inhalation [31]. The corresponding pressure was applied at the outlet of the fluid domain [32] (see Figure 4). Properties of air were assumed to be those at 27°C. The airways were assumed to be of a homogeneous material with a density of 1365.6 kg/m3 [33], a Young’s modulus of elasticity in longitudinal direction of 130.89 kPa [34], a Young’s modulus of elasticity in circumferential direction of 74.07 kPa [35], and Poisson’s ratio of 0.45 [35]. A zero-displacement boundary condition was applied at the airways on both inlet and outlet of the solid domain to represent a tethering of the airways from other tissues and organs [36]. A no-slip boundary condition was defined at the fluid-structure interface (see Figure 3).

2.6. Methods of Analysis

Airflow velocity, airway pressure, and airway stresses were calculated in this study. Stresses in the airways during inhalation are from fluid shear force and airway pressure. Fluid shear force during inhalation creates wall shear stress (WSS) at the surface of the airways. WSS is the tangential stress at a wall due to fluid viscosity and is related to a transverse velocity gradient [37]. In contrast, airway pressure acting in the normal direction to the airways creates stresses across thickness of the airways. Stresses from the airway pressure are normal and shear stress distributions in all directions. In this study, stress distributions in the airways were analyzed employing the longitudinal and von Mises stresses. The von Mises stress is an average stress in all directions and is associated with distorting the shape of material. Any material will yield if the von Mises stress is greater than its yield strength [38].

Effects of airway tissue flexibility on airflow velocity, airway pressure, and WSS in the bifurcation were analyzed assuming the airways to be rigid or flexible. For the rigid model, analysis was performed only on the fluid domain. The bifurcation in this case acts like a rigid tube that cannot be deformed by fluid forces from the airflow; therefore, there is no stress in the airways. For the flexible model, analysis was performed on both fluid and solid domains. The bifurcation in this case acts like a flexible tube that can be deformed by fluid forces from the airflow. Therefore, stresses in the airways are considered in this analysis. The FSI analysis [20] was implemented for the flexible model.

To investigate the effect of tissue flexibility on stress distributions in the airways, simulations were performed using orthotropic and isotropic material models. For the orthotropic material model, the properties of the airways described in Section 2.5 were used. In contrast to the orthotropic material model, the Young’s modulus of elasticity of 130.89 kPa [34] and Poisson’s ratio of 0.45 [35] was used for the isotropic material model.

3. Results and Discussion

The distributions of airflow velocity, airway pressure, WSS, longitudinal stress, and von Mises stress for the flexible model were evaluated only with orthotropic material model as the distributions for both flexible and rigid models are similar. The effect of tissue flexibility on airflow velocity, airway pressure, and WSS and the effect of material models on longitudinal and von Mises stresses are also discussed.

3.1. Airflow Velocity

The distributions of airflow velocity in the airways during the peak airflow are shown in Figure 5. The maximum airflow velocity during the peak airflow was 3.692 m/s. High airflow velocity spread throughout G3. Airflow velocity in the medial side of G4 was higher than that in the lateral side of G4 because the secondary flow after the bifurcation moved air toward the medial side of G4. Airflow velocity in both G3 and G4 was symmetric; however velocity profiles in G5 were not symmetric. Airflow velocity in branch G5M was higher than that in branch G5L because of the high airflow velocity in the medial side of G4.

3.2. Airway Pressure and von Mises Stress

The distributions of airway pressure and von Mises stress in the airways during the peak airflow are shown in Figure 6. The maximum airway pressure during the peak airflow was −3.209 Pa. High airway pressure was observed at the beginning of G3 and the bifurcations. The pressure at the first bifurcation between G3 and G4 was about two times higher than that at the second bifurcation between G4 and G5 since the airflow that impinged the first bifurcation had more momentum (higher velocity) than the air that impinged the second bifurcation. The von Mises stress in the airways, which resulted from the airway pressure, was observed to be high at the beginning of G3 and the bifurcations. The von Mises stresses in branch G5M were higher than those in branch G5L. The maximum von Mises stress during the peak airflow was 1142 Pa.

3.3. Wall Shear Stress (WSS) and Longitudinal Stress

The distributions of wall shear stress (WSS) at the airways and longitudinal stress in the airways during the peak airflow are shown in Figure 7. The maximum WSS, 0.4185 Pa, occurred at the bifurcations since the velocity gradient is highest at those locations. WSS in the medial side of G4 was higher than that in the lateral side of G4, and WSS in G5M was higher than that in G5L. When the locations of the high WSS are compared with the particle deposition locations from the previous study by Longest and Vinchurkar [26], we can see that the particles deposit in the high WSS locations (see Figure 8). The longitudinal stress in the airways, which resulted from the airway pressure, was observed to be high where WSS was observed to be high. The maximum longitudinal stress during the peak airflow was −1497 Pa. The negative value of the longitudinal stress indicates that the length of the airways decreased during the peak airflow.

3.4. Tissue Flexibility Effect

Figure 9 shows the effect of tissue flexibility on the airflow velocity, airway pressure, and WSS during inhalation. As can be seen from this figure, tissue flexibility affected each flow characteristic to a different degree. Airway pressure and WSS were significantly affected by tissue flexibility. The maximum differences between the flexible and rigid models were 2%, 7%, and 6%, for airflow velocity, airway pressure, and WSS, respectively. The airflow velocity from the flexible model was higher than that from the rigid model for the first 0.9 seconds since the airways contracted due to the negative pressure inside the airways. However, the airflow velocity from the flexible model was lower than that from the rigid model for the last 0.9 seconds since the airways expanded due to the positive pressure inside the airways. The airway pressure from the flexible model was lower than that from the rigid model throughout the inhalation process. The movement of the airway tissue decreased the airway pressure, and the maximum differences between the flexible and rigid models were highest at the peak airflow. The increase in airway pressure associated with the rigid wall suggests that people with stiff airways, for example, elderly people [39] or asthma patients [40] can experience high airway pressure during normal inhalation. Tissue flexibility decreased WSS for the first 0.9 seconds; however, it increased WSS for the last 0.9 seconds. The effects of tissue flexibility on WSS from this study were similar to results from the previous studies in the abdominal aortic aneurysm (AAA) by Leung et al. [41] and Torii et al. [42] as well as Scotti and Finol [43]. Their results showed that tissue flexibility can increase or decrease WSS, and the influence on WSS is highly depended on AAA geometry.

3.5. Material Model Effect

Figure 10 shows the effect of the material models on the longitudinal and von Mises stresses. The longitudinal stress from the orthotropic material model was lower than that from the isotropic material model for entire inhalation since the isotropic material model was stiffer than the orthotropic material model. The differences between both material models for the longitudinal stress are in the ranges of 22–52%. The differences between both material models were observed to be lowest near the peak airflow. The von Mises stress from the orthotropic material model was higher than that from the isotropic material model for the first 0.9 seconds; however, the von Mises stress from the orthotropic material model was lower than that from the isotropic material model for the last 0.9 seconds. The differences between both material models for the von Mises stress are in the ranges of 25–28%. The maximum difference was observed near the peak airflow.

4. Limitations

In this study, the airway geometry was based on an idealized ICRP [22] symmetric model. However, a study by Horsfield et al. [44] showed that airway diameters and branching airways were asymmetric. In addition, material properties of the airways in the present study were assumed to be linear. Studies by Ito et al. [45] and Smith et al. [46] showed that the airways exhibited viscoelastic properties and nonlinear dynamic behaviors. Further study is needed to investigate the effect of the airway diameters, branching angles, and airway properties on flow characteristics and stresses in the airways.

5. Conclusions

The airflow velocity, airway pressure, WSS, and stresses within the airway generations 3 to 5 were analyzed in this study using the finite-element method with the FSI algorithm. The analysis was performed to investigate the effects of tissue flexibility and the material model on flow characteristics and stresses in the airways during inhalation. The simulation results showed that tissue flexibility decreased the airway pressure and altered airflow velocity and WSS. The simulation results also showed that the material model of the airways significantly affected stresses in the airways. The results from this study highlight the importance of incorporating tissue flexibility along with orthotropic properties into a computational model that is developed to study flow characteristics and stresses in the airways.