Introduction

Numerical studies in the field of the respiratory system provide useful information about the respiratory flow structures, deposition patterns of particulate matter, and pharmaceutical aerosols. Since the beginning of the COVID-19 pandemic, studies in this field have become more important for researchers to assess the behavior of viruses entering the respiratory tracts [ 1 ]. In Computational Fluid Dynamics (CFD), the design of geometry and grid generation is considered preprocessing stages with an important role in the achievements of numerical research. Numerical studies in the field of human airways are accompanied by limitations, such as computational cost and limited resolution of diagnostic medical imaging equipment [ 2 ]. As a result, respiratory flow in truncated geometries should be always evaluated.

Different methods are considered for producing the airways model, which is associated with limitations. The airways model can be produced in three main methods: based on idealized models [ 3 ], based on stochastic methods [ 4 ], and based on medical imaging equipment [ 5 , 6 ]. Researchers have compared the respiratory flow in realistic and idealized models from different aspects or have used idealized models to study a particular parameter. Heidarinejad et al. studied the effect of the cartilaginous tracheal rings on the airflow structures by modifying the Horsfield geometry [ 7 ], and also observed small recirculating structures between ring cavities, leading to significant increases in the wall shear stress. Since the wall shear stress reaches its peaks at cartilaginous rings, the maximum chances of wall injuries, such as epithelial damage are in these areas [ 8 ]. Therefore, wall curvature in realistic models has a significant effect on respiratory flow structures and identifies vulnerable areas [ 9 ]. Johari et al. compared the velocity and pressure data in realistic, simplified, and oversimplified models in terms of roughness and curvature on the wall [ 10 ] and found significant deviations at high flow rates in regions with high geometric complexities, especially in areas with rapid cross-sectional changes. However, idealized and stochastic models can provide a complete geometry of the entire lung [ 11 ], the flow characteristics obtained from these models do not fit well with the results of the realistic models [ 12 ]. Furthermore, a high-computation cost is imposed to generate a high-resolution grid due to the small size of the last generations of the lungs [ 13 ]. Therefore, optimal geometry is considered an important challenge for researchers. Wall roughness and small vortex structures near the wall affect the particle deposition pattern in the human upper airways [ 14 ]. Hence, realistic geometries lead to more accurate information about the respiratory flow structures and particle deposition patterns. On the other hand, specifications of the middle and last generations of the lungs are not accessible by diagnostic medical imaging equipment [ 15 ]. As a result, a truncated geometry of the lungs is always investigated.

A limited number of studies investigated the effect of geometry truncation on flow variables. Qi et al. compared the spatial and temporal variables of the flow field based on two models with 122 and 22 outlets [ 16 ] and also showed that the lobar air distribution considerably correlated with the outlet area ratio. Moreover, the upper respiratory tract was not considered with a significant effect on the flow characteristics. Choi et al. evaluated the effect of the upper airway truncation in four different levels [ 17 ] and also developed a proper inlet boundary condition for the truncated sections. According to their results, the respiratory flow around the larynx is highly affected by the laryngeal jet. Lin et al. and Xi et al. showed that the upper airway truncation, such as oropharynx and larynx, eliminated velocity parabola and turbulence in the trachea, leading to different flow structures [ 5 , 9 ]. Additionally, Xiao et al. concluded that subglottic inflow truncation delayed the location of the jet breakup issuing from the tracheal contraction [ 18 ]. Tena et al. developed a boundary condition for the truncated sections by applying velocities from similar locations of an extended branch [ 2 ]. This method can only be used for idealized geometries, and an extended branch’s fully resolved airflow field must be available. Since fluid properties diffuse in all directions, downstream areas of the computational domain always affect upstream areas [ 19 ]. Hariprasad et al. showed that flow patterns in the central airway were highly affected by obstructions in the lower airways [ 20 ]. The present study aimed to evaluate the effect of bronchial airways, as a major resistance to the airflow path, on the respiratory flow characteristics of upstream areas. For this purpose, geometries with zero-generation and 5^th-generation truncation levels were employed to assess the impact of bronchial airways truncation on tracheal airflow characteristics.

Material and Methods

1. Airway Modelling and Description of Mesh Parameters

In this numerical study, Computed Tomography (CT) images were employed to produce a realistic model of the human respiratory airways belonging to a 35-year-old healthy man. The Digital Imaging and Communications in Medicine (DICOM) file contains 485 images with a spatial resolution of 0.6×0.6 mm², including the upper respiratory tract (nasal cavities, pharynx, and larynx) and the lower airways until the 5^th generation (Figure 1). The airflow characteristics in a truncated model of the original geometry were also investigated. The truncated model is up to the zero generation (up to the end of the trachea).

Figure 1. The realistic human airways model from nostrils up to the 5^th generation

The airways models were meshed by unstructured tetrahedral elements with a maximum cell size of 0.5 mm using ANSYS ICEM-CFD. The near-wall grid size must sufficiently be refined as a high-velocity gradient occurs near the wall, especially in the viscous sub-layer [ 21 ]. Six prism layers were attached to the airways wall with the first cell height of 0.05 mm and a growth rate of 1.1, as shown in Figure 1. Consequently, the first layer closest to the wall was located in the Y-Plus range of less than one in almost all areas of the computational domain. Accordingly, the computational grid well captures all fluidic layers near the wall, including the viscous sub-layer, the buffer layer, and the logarithmic layer. As a result, the correct behavior of the wall returns to the computational domain. A high-quality mesh results in an accurate and stable numerical simulation, especially in turbulent wall-bounded flows. A total of 4.2 and 6.3 million computational cells were generated for geometries with zero-generation and 5^th-generation truncation levels, respectively.

2. Governing Equations

The computational fluid dynamics techniques was employed to solve the airflow governing equations. A finite volume solver (ANSYS Fluent) was employed to solve the unsteady, incompressible, and three-dimensional Navier-Stokes equations using the large eddy simulation (LES) turbulent technique coupling with the WALE sub-grid scale model. The large eddy simulation technique is validated for solving transitional to turbulent flows, [ 22 ] which has a lower computational cost than the direct numerical simulation (DNS) method [ 23 ]. LES has an accurate performance in predicting separated flows with high levels of unsteadiness [ 24 ]. Furthermore, the WALE sub-grid scale model can capture the transition process and return a zero turbulent viscosity in laminar shear flows [ 25 ]. In the WALE model, eddy viscosity is determined without explicit filtering and instead employs a tensor invariant to apply the strain and rotation effects. Additionally, this model is applicable to wall-bounded flows and complex geometries due to the use of only local gradients to compute the eddy viscosity and return the correct wall asymptotic behavior [ 26 ].

The filtered continuity and momentum equations for unsteady, incompressible, and isothermal airflow fields are expressed as follows:

$\frac{\partial \stackrel{-}{U_i}}{\partial x_i}=0$ (1)

$\frac{\partial \stackrel{-}{U_i}}{\mathrm{\partial t}}+\frac{\partial }{\partial x_j}\left(\stackrel{-}{U_i}\stackrel{-}{U_j}\right)=-\frac{1}{\rho }\frac{\partial \stackrel{-}{P}}{\partial x_i}+\frac{\partial }{\partial x_j}\left[\upsilon \right(\frac{\partial \stackrel{-}{U_i}}{\partial x_i}+\frac{\partial \stackrel{-}{U_j}}{\partial x_i}\left)\right]-\frac{\partial {\tau }_{\mathrm{ij}}}{\partial x_j}$ (2)

The over-bar is the symbol of filtering, and $\stackrel{-}{U_j}$, $x_i$, and $\stackrel{-}{P}$ are the filtered velocity component, position component, and filtered static pressure, respectively. ρ is the fluid density assumed 1.2 kg.m^-2. ${\tau }_{\mathrm{ij}}$ is the sub-grid scale stresses tensor related to the large-scale strain rate tensor through the following Equations:

$({\tau }_{\mathrm{ij}}-\frac{{\delta }_{\mathrm{ij}}}{3}{\tau }_{\mathrm{kk}})=-2{\upsilon }_{\mathrm{SGS}}\stackrel{-}{S_{\mathrm{ij}}}$ (3)

$\stackrel{-}{S_{\mathrm{ij}}}=\frac{1}{2}(\frac{\partial \stackrel{-}{U_i}}{\partial x_j}+\frac{\partial \stackrel{-}{U_j}}{\partial x_i})$ (4)

The sub-grid scale viscosity ${\upsilon }_{\mathrm{SGS}}$ is modeled as [ 25 ]:

${\upsilon }_{\mathrm{SGS}}={(C_w\Delta )}^2\frac{{\left(S_{\mathrm{ij}}^d{S}_{\mathrm{ij}}^d\right)}^{3/2}}{{\left(\stackrel{-}{S_{\mathrm{ij}}}\stackrel{-}{S_{\mathrm{ij}}}\right)}^{5/2}+{\left(S_{\mathrm{ij}}^d{S}_{\mathrm{ij}}^d\right)}^{5/4}}$ (5)

$s_{\mathrm{ij}}^d=\frac{1}{2}(\stackrel{-}{g_{\mathrm{ij}}^2}+\stackrel{-}{g_{\mathrm{ji}}^2})-\frac{1}{3}{\delta }_{\mathrm{ij}}\stackrel{-}{g_{\mathrm{kk}}^2}$ (6)

Δ is the LES sub-grid length scale $(\Delta ={\left(\mathrm{\Delta X\Delta Y\Delta Z}\right)}^{1/3})$, and $s_{\mathrm{ij}}^d$ is the traceless symmetric part of the square of the velocity gradient tensor; $\stackrel{-}{g_{\mathrm{ij}}}$ is velocity gradient tensor $(\stackrel{-}{g_{\mathrm{ij}}}=\frac{\partial \stackrel{-}{U_i}}{\partial x_j})$; the WALE constant $C_w$ is considered 0.325.

3. Numerical Setup

The gradient terms were approximated by the least-squares cell-based differencing scheme, which is a highly accurate method in skewed networks with a lower computational cost than other methods. Convection and transient terms were discretized by bounded central differencing and bounded second-order implicit differencing methods, respectively, in order to prevent increasing the numerical diffusion error. The time step size is set to 10^-4 s, which ensures the Courant-Friedrichs-Lewy (CFL) number of less than 0.8 throughout the computational domain and improves the stability as well as the accuracy of the numerical simulation.

The convergence of the averaged velocity field was evaluated as a criterion to access the statically stable solution. The time interval of 0.1 s is considered for averaging since instantaneous effects of the airflow field disappear in this averaging process. Numerical simulations converged in the geometry after 0.5 and 0.4 s with the 5^th-generation truncation level and the zero-generation truncation level, respectively.

4. Boundary Conditions

4.1. Inlets Boundary Condition

Different methods are used for defining the stimulus boundary condition in numerical simulations, such as using experimental data at the boundaries [ 27 ] or the stimulus boundary condition at the end of the airways model [ 28 ] since the respiration work is done by the movement of the diaphragm vertically at the end of the lungs. Applying these methods to realistic geometries with bronchial airways is accompanied by some limitations.

In this paper, the uniform velocity profile is employed at the nose inlets to provide an equal situation for comparing the results of the two geometries. The velocity magnitude at the inlets is adjusted according to the flow rate of 45 L/m, showing the heavy breathing condition; the inlet Reynolds number is about 4200. The inlet velocity profile affects the flow pattern in idealized models [ 29 ]; however, the effect of the inlet velocity profile disappears at a short distance from the nose inlets in realistic models due to severe cross-sectional changes in the nasal cavity [ 15 ]. Therefore, the uniform velocity assumption has no effect on the airflow characteristics in the trachea as the intended area.

4.2. Outlets Boundary condition

One of the major challenges in this field of study is defining an appropriate outlet boundary condition, especially in geometries with branching structures. The branching structure of the lungs affects the pressure in the entire geometry and acts as a resistance to the airflow path. Alzahrany & Banerjee developed a pressure-based outlet boundary condition [ 30 ], which was a function of airways resistance. The procedure of calculations is based on the Weibel geometry, which is an idealized model. Comerford et al. proposed a structured tree impedance outflow boundary condition to simulate the airflow in patient-specific human lungs [ 31 ]. The outlet pressure value varies based on healthy or diseased lungs, age, and respiration states. Moreover, the uniqueness of image-based geometries complicates the situation. On the other hand, the pressure value varies in different branches of the same generation.

In this study, the pressure outlet and outflow boundary conditions were compared, in which the pressure outlet boundary condition showed a better adaptation to the physics of respiration. The pressure outlet boundary condition employed zero static (gauge) pressure. This method extrapolated the interior cells’ stress and all other flow quantities.

In the outflow boundary condition, a zero-diffusion flux was applied for all flow variables; the flow rate weighting obtained from the pressure outlet boundary condition was used for the outflow boundary condition.

4.3. Wall Boundary condition

All velocity components at the wall are set to zero in order to apply the no-slip boundary condition.

5. Validation

The numerical simulation procedure has been validated by Jayaraju et al.’s experimental and numerical results [ 32 ]. They also developed a simplified model of the upper respiratory tract considering this region’s curvatures, expansions, and contractions, as shown in Figure 2. The numerical simulation was performed using the LES approach coupling with the WALE sub-grid scale model with 3.2 million computational cells. The velocity magnitude at the inlet is adjusted according to the Reynolds number of 2527, which is equal to the value considered in the experiment.

Figure 2. a) section A, b) section B – PSD (Power Spectral Density) analysis at two points in the center, c) section 3, d) section 5, and e) The LES (Large Eddy Simulation) index of quality at different sections

The time-averaged velocity data of the present simulation were employed for comparison with the experimental and numerical data available from Jayaraju et al. Figure 2 (a, b) shows the proper adaptation of the mean velocity profiles of the present study with the particle image velocimetry (PIV) measurements and numerical data of Jayaraju et al. at specified sections. Therefore, the numerical simulation procedures used in this section, such as the properties of the generated mesh, time step size, and discretization schemes, were also used to study the airflow characteristics in the proposed realistic geometry.

Results

1. Numerical Simulation Accuracy

The simulation of turbulent flows over complex geometries using unsteady approaches, such as DNS and LES, is accompanied by complexities originating from turbulence’s nature [ 33 ]. The energy spectrum law is a criterion for investigating the accuracy of turbulent flows simulation, derived theoretically by Kolmogorov [ 34 ], who illustrated that turbulent flows are locally homogeneous and isotropic at sufficiently high Reynolds numbers. Moreover, based on Kolmogorov’s results, turbulent flows must be statistically equilibrium in the range of high wavenumbers [ 35 ]. The energy spectrum law (Kolmogorov’s K^-5/3 law) is defined as:

$E\left(K\right)=K_0{\epsilon }^{2/3}{K}^{-5/3}{\lambda }^{-1}\ge K\ge l^{-1}$ (7)

where $E\left(K\right)$ is the scalar energy distribution function per mass unit for a certain wave number value K (the Fourier transform of the kinetic energy); $K_0$ is the Kolmogorov number, which is between 1.4 and 1.8. Additionally, ε, λ, and l are the dissipation rate, the characteristic length scale for the viscous eddies, and the integral length scale, respectively.

Figure 2 (c, d) shows that the energy spectrum density matches the K^-5/3 slope at high wavenumbers (Taylor scale) using the frequency analysis of instantaneous velocity data, i.e., the transmission of energy from large-scale structures to small-scale structures (Kolmogorov scale) was accurately modeled, highly dependent on the number of computational cells and its quality.

The LES index of quality (LES-IQ) also is a criterion, measuring the resolved turbulent kinetic energy ratio to the total turbulent kinetic energy. In an accurate numerical simulation using LES models, 80% of the turbulent kinetic energy should be determined [ 7 ]. The LES_IQ_ν relates the turbulent viscosity to the laminar viscosity through the following equation [ 36 ]:

${\mathrm{LES\_IQ}}_{\upsilon }=\frac{1}{1+0.05{\left(\frac{⟨{\upsilon }_{\mathrm{t.eff}}⟩}{\upsilon }\right)}^{0.53}}$ (8)

$⟨{\upsilon }_{\mathrm{t.eff}}⟩$ and υ are the time-averaged effective turbulent viscosity and the laminar viscosity (fluid viscosity).

According to Figure 2 (e), the minimum value of the LES_IQ_ν in different sections of the computational domain is greater than 0.9, showing the amount of the resolved turbulent kinetic energy is in an acceptable range. The PSD analysis and LES-IQ results were obtained from the numerical simulation performed in the geometry with the 5^th-generation truncation level. Since the results of these two criteria indicate the high accuracy of the numerical solution and the appropriate quality and quantity of the generated mesh, the numerical simulation procedure performed in this geometry was also used for the geometry with the zero-generation truncation level.

2. Airflow Characteristics

According to the physics of the respiratory flow, the branching structure of the lungs can act as a resistance to the flow path and increase the pressure throughout the geometry. Considering the first few generations of the lungs will prevent abnormal pressure drops in the trachea or upstream areas. Moreover, the outlet boundary condition typically affects the pressure drop throughout the geometry.

Figure 3 (a) illustrates the effect of the pressure outlet and outflow boundary conditions on the pressure drop in geometries with zero-generation and 5^th-generation truncation levels. As shown in Figure 3 (a), the pressure drop in the two mentioned geometries using the outflow boundary condition is approximately the same, while a 20-to-25-unit difference is observed using the pressure outlet boundary condition. Therefore, the pressure outlet boundary condition provides a situation that is more consistent with the physics of respiration. In addition, a logical pressure drop is observed across the geometry using the pressure outlet boundary condition compared to experimental data [ 37 ]. According to the mentioned points, the pressure outlet boundary condition was employed to investigate the effect of bronchial airways truncation on tracheal airflow characteristics in the rest of this paper.

Figure 3. Mean pressure drop using pressure outlet and outflow boundary conditions, b) normalized turbulent kinetic energy (TKE/ U²), and c) vorticity Magnitude in geometries with zero-generation and 5^th-generation truncation levels

Evaluating flow instabilities can lead to examining the impact of bronchial airway truncation. For the first step, normalized turbulent kinetic energy was investigated in different slices in Figure 3 (b). Section 8 in Figure 3 (b) (Longitudinal section of the trachea) shows that the airflow field is more turbulent in the geometry with the zero-generation truncation level, especially in the distance from the larynx to the end of the trachea. Detailed analysis of the transverse sections showed that the turbulent kinetic energy around the pharynx grows abnormally (Figure 3 (b, section 3)).

Another important parameter analyzed is vorticity, which is a microscopic parameter representing the local rotational motion in the computational domain. Figure 3 (c) shows the vorticity magnitude at different levels of the two geometries. The difference in vorticity magnitude between geometries with and without branching structure has begun to grow from section 3 (around the pharynx) and reaches its peak in sections 4, 5, and 6, as shown in Figure 3 (c). In the mentioned sections, the vorticity magnitude is significantly higher in the geometry with the zero-generation truncation level.

Flow patterns using tracheal streamlines are shown in Figure 4 (a, b); the colors of streamlines indicate the Z and U components of the velocity. Since the direction of primary flow in the trachea is in the negative direction of the Z-axis, secondary flows and vortices are obviously detectable by the red color.

Figure 4. Tracheal streamlines colored: a) the Z component, b) the U component of the velocity, c) sectional streamlines at different sections of the trachea. Forward streamlines colored by the Z component of the velocity around, d) section 4, and e) section 5

Figure 4 (a, b) shows that more robust secondary flows are observed above and below the larynx in the 5^th-generation truncation level geometry.

Forward streamlines in sections 4 and 5 are shown in Figure 4 (d, e) to evaluate the airflow structure around the larynx in detail. In Figure 4 (d, around section 4), a more extended vortex was observed in the geometry with the 5^th-generation truncation level. In Figure 4 (e, around section 5), the vortex formed in the lower area of the larynx is broader and longer in the geometry with the 5^th-generation truncation level.

In Figure 4 (e, around the section 5), the vortex formed in the lower area of the larynx is wider and longer in the geometry with the 5^th-generation truncation level.

In Figure 4 (c), streamlines are shown at different sections between the larynx and the end of the trachea. Streamlines in all three sections represent the existence of wider vortex structures in the geometry with the 5^th-generation truncation level.

The structure of vortices affects the wall shear stress locally, which is generally more severe at wall shear stress peaks. Figure 5 shows the time-averaged wall shear stress at different sections along the trachea.

Figure 5. The variation of wall shear stress using geometries with zero generation and 5^th-generation truncation level

According to Figure 5, the wall shear stress has reached higher peaks in the geometry with the 5^th-generation truncation level. A 14-to-31% difference is observed between the major peaks of the wall shear stress using the mentioned geometries.

Discussion

This paper aimed to investigate the effect of bronchial airways truncation on the airflow characteristic as an important factor in the accuracy of numerical simulations in this field of study. For this purpose, the respiratory flow was studied in geometries with and without bronchial airways.

The turbulent kinetic energy has reached higher peaks in the geometry with the zero-generation truncation level, as shown in Figure 3 (b). The absence of the branching structure of the lungs leads the air to flow more freely in the geometry with the zero-generation truncation level. Therefore, factors producing perturbation, such as tracheal wall curvatures, create more severe fluctuations in the velocity field of this geometry, resulting in the instantaneous velocity field deviating more intensely from the mean velocity field. The turbulent kinetic energy in the geometry, including the branching structure of the lungs is intensified at a shorter distance after the larynx (Figure 3 (b, section 8-L1)), while this distance is slightly greater in the geometry with the zero-generation truncation level (Figure 3 (b, section 8-L2)). Studies that employed geometries with bronchial structures have confirmed the intensification of the turbulent kinetic energy immediately after the larynx [ 38 - 40 ]. There is not any noticeable difference in the intensity of turbulent kinetic energy around the nasal cavity and the pharynx between the two geometries, according to the turbulent kinetic energy (TKE) in sections 2 and 7 in Figure 3 (b). The airflow characteristics in the nasal cavity and adjacent areas are strongly influenced by the high-curvature structure of the nasal geometry. Moreover, entering the airflow into the bronchi in the extended model has increased the turbulent kinetic energy just at the end of the trachea and creates a relative balance in the intensity of turbulent kinetic energy in this area between geometries with and without bronchial airways. Furthermore, according to Figure 3 (b, section 7), the core of the turbulent kinetic energy is transferred to the center of the section in the geometry with the zero-generation truncation level.

Vorticity magnitude has also increased along the trachea. As mentioned above, the absence of the branching structure of the lungs has led the air to flow more freely in the trachea, and factors producing perturbation create more severe fluctuations in the velocity field that these fluctuations in the airflow field have caused more severe slipping of fluidic layers on each other, which ultimately reinforces the tendency of each point of the computational domain to rotate. As the airflow exhibits more chaotic behavior in the geometry with the zero-generation truncation level, the probability of forming larger-scale vortices is reduced. The vorticity has a completely different condition at the end of the trachea. Unlike upper areas of the trachea, the vorticity is reinforced adjacent to the bronchi in the geometry with the 5^th-generation truncation level (Figure 3 (c, section 7)). Changing the direction of airflow from the trachea to the bronchi and cross-sectional changes in this region are influential factors in the vorticity, resulting in more severe fluctuations in the velocity field of this region.

The upper and lower areas of the larynx are associated with significant increases in the cross-sectional area compared to the larynx, and act as a cavity to the flow path; large-scale vortices form in these parts. A comparison of vortex properties around the larynx to the end of the trachea can obviously show the effect of bronchial airway truncation. According to Figure 4 (b-L2), a stronger and larger-scale vortex is generated in the geometry with the 5^th-truncation level. The differences in vortex structures are much more severe in the lower part of the larynx. The vortex formed in the geometry with the 5^th-generation truncation level is much larger than the vortex of the same location in the geometry with the zero-generation truncation level (Figure 4 (a-L1)). Larger-scale vortices have also formed at the end of the trachea in the geometry with the 5^th-generation truncation level. Drastic changes in the direction of airflow entering the bronchial airways can be an effective factor in the formation of larger-scale vortices at the end of the trachea in the geometry, including the branching structures of the lungs.

These larger-scale vortices in the extended model have caused the primary flow passing next to the vortices to accelerate more intensely, leading to increases in the local velocity field eventually. This issue can be clearly observed by the color of the primary flow passing next to the displayed vortices in Figure 4 (d, e). The dimensions of vortices directly affect the local wall shear stress rate, which is an important parameter in determining vulnerable areas. Therefore, accurate calculation of this parameter is vital.

The right and left main bronchi greatly influence the formation of secondary flows at the end of the trachea and at the entrance of the main bronchi.

According to Figure 1, the right bronchus of the lung has a steeper angle to the trachea than the left bronchus, causing a larger-scale vortex to form at the entrance of the right bronchus and at the end of the trachea, which is the result of the airflow separation in this area (Figure 4 (c, section 7)).

The vortex structure affects the particle deposition pattern and its residence time. Particle residence time increases during passing from larger-scale vortices, and the probability of particle redirection increases. Since the airflow in respiratory airways is severely bounded to the wall, the lack of accurate simulation would result in changes in determining the particle deposition pattern.

A significant number of numerical studies in the field of human respiratory tracts investigated airway obstructions and their effects on the airflow characteristics [ 20 ]. Wall shear stress is a key determinant in identifying vulnerable areas. Epithelium cells located at the inner side of the airway wall are sensitive to the wall shear stress rate. As a result, accurate calculation of wall shear stress can be very effective in identifying vulnerable areas and estimating the extent of the damage. In this study, significant differences are observed in the shear stress values in the two geometries with the zero-generation truncation level and the extended one. According to Figure 5, the major peak of each chart is located adjacent to areas of the computational domain with high local velocity magnitude. Each vortex obstructs the airflow path, resulting in the primary flow passing next to the vortex to accelerate and eventually increasing the local velocity field. Larger-scale vortices result in the primary flow accelerating more intensely. Since larger-scale vortices formed in the geometry with the 5^th-generation truncation level, the wall shear stress related to this geometry has reached higher peaks.

In this study, limitations, such as high computational cost and limited hardware facilities caused to prevent studying the airflow characteristics during a complete respiratory cycle. In addition, the limited spatial accuracy of medical imaging equipment made it impossible to study a complete geometry.

Conclusion

In this paper, the effect of bronchial airways truncation was assessed on tracheal airflow characteristics using geometries with zero generation and 5^th-generation truncation levels. Some of the most important achievements of this paper are as follows: 1) Higher level of turbulent kinetic energy is observed along the trachea in the geometry with the zero-generation truncation level. The elimination of bronchial airways of the lungs leads the air to flow more freely in this geometry. Therefore, factors producing perturbation, such as tracheal wall curvatures, create more severe fluctuations in the velocity field, 2) the more intense fluctuations in the geometry without branching structure increased the tendency of each point to rotate and reinforced the vorticity along the trachea in the geometry with the zero generation truncation level, and 3) since the airflow exhibits a more chaotic behavior in the geometry with the zero-generation truncation level, larger-scale vortices are less likely to form. On the other hand, larger-scale vortices have formed in the geometry with the 5^th-generation truncation level, and 4) the primary flow passing next to the larger-scale vortices accelerates more intensely, resulting in increases in the local velocity field. As a result, the wall shear stress has reached higher peaks in the geometry with the 5^th-generation truncation level.

Authors’ Contribution

S. Farhoodi and Dr. MH. Roozbahani conducted the numerical simulations and analyzed the results of this work. S. Farhoodi wrote the original draft of this paper. The project was supervised by Professor Gh. Heidarinejad. All the authors have read and approved the final manuscript of this article.

Ethical Approval

The Iranian Research Institute for Information Science and Technology (IranDoc) approved the protocol of this research study (Code: 2665014).

Informed consent

The CT images have been provided to the researcher with informed consent.

Conflict of Interest

None

References

1. Vuorinen V, Aarnio M, Alava M, Alopaeus V, Atanasova N, et al. Modelling aerosol transport and virus exposure with numerical simulations in relation to SARS-CoV-2 transmission by inhalation indoors. Saf Sci. 2020; 130:104866. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
2. Tena AF, Fernández J, Álvarez E, Casan P, Walters DK. Design of a numerical model of lung by means of a special boundary condition in the truncated branches. Int J Numer Method Biomed Eng. 2017; 33(6):e2830. DOI | PubMed
3. Horsfield K, Cumming G. Angles of branching and diameters of branches in the human bronchial tree. Bull Math Biophys. 1967; 29(2):245-59. DOI | PubMed
4. Koblinger L, Hofmann W. Stochastic morphological model of the rat lung. Anat Rec. 1988; 221(1):533-9. DOI | PubMed
5. Lin CL, Tawhai MH, McLennan G, Hoffman EA. Characteristics of the turbulent laryngeal jet and its effect on airflow in the human intra-thoracic airways. Respir Physiol Neurobiol. 2007; 157(2-3):295-309. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
6. Roozbahani MH, Heidarinejad M, Heidarinejad G. Impacts of human tracheal cartilaginous rings on tracheobronchial flow structures. 15th Conf Int Soc Indoor Air Qual Clim INDOOR AIR; Philadelphia, Pennsylvania, USA: Drexel University; 2018. p. 1-4.
7. Heidarinejad G, Roozbahani MH, Heidarinejad M. Studying airflow structures in periodic cylindrical hills of human tracheal cartilaginous rings. Respir Physiol Neurobiol. 2019; 266:103-14. DOI | PubMed
8. Srivastav VK, Paul AR, Jain A. Effects of cartilaginous rings on airflow and particle transport through simplified and realistic models of human upper respiratory tracts. Acta Mech Sin. 2013; 29(6):883-92. DOI
9. Xi J, Longest PW, Martonen TB. Effects of the laryngeal jet on nano- and microparticle transport and deposition in an approximate model of the upper tracheobronchial airways. J Appl Physiol. 2008; 104(6):1761-77. DOI | PubMed
10. Johari NH, Osman K, Helmi NH, Abdul Kadir MA. Comparative analysis of realistic CT-scan and simplified human airway models in airflow simulation. Comput Methods Biomech Biomed Engin. 2015; 18(1):48-56. DOI | PubMed
11. Weibel ER. Geometric and dimensional airway models of conductive, transitory and respiratory zones of the human lung. In: Morphometry of the human lung; Berlin, Heidelberg: Springer; 1963.
12. Bocanegra Evans H, Castillo L. Index-matched measurements of the effect of cartilaginous rings on tracheobronchial flow. J Biomech. 2016; 49(9):1601-6. DOI | PubMed
13. Koullapis PG, Hofemeier P, Sznitman J, Kassinos SC. An efficient computational fluid-particle dynamics method to predict deposition in a simplified approximation of the deep lung. Eur J Pharm Sci. 2018; 113:132-44. DOI | PubMed
14. Shi H, Kleinstreuer C, Zhang Z. Modeling of inertial particle transport and deposition in human nasal cavities with wall roughness. Journal of Aerosol Science. 2007; 38(4):398-419. DOI
15. Koullapis PG, Kassinos SC, Bivolarova MP, Melikov AK. Particle deposition in a realistic geometry of the human conducting airways: Effects of inlet velocity profile, inhalation flowrate and electrostatic charge. J Biomech. 2016; 49(11):2201-12. DOI | PubMed
16. Qi S, Zhang B, Teng Y, Li J, Yue Y, Kang Y, Qian W. Transient Dynamics Simulation of Airflow in a CT-Scanned Human Airway Tree: More or Fewer Terminal Bronchi? Comput Math Methods Med. 2017; 2017:1969023. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
17. Choi J, Tawhai MH, Hoffman EA, Lin CL. On intra- and intersubject variabilities of airflow in the human lungs. Phys Fluids. 2009; 21(10):101901. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
18. Xiao Q, Cetto R, Doorly DJ, Bates AJ, et al. Assessing Changes in Airflow and Energy Loss in a Progressive Tracheal Compression Before and After Surgical Correction. Ann Biomed Eng. 2020; 48(2):822-33. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
19. Crowe CT. Multiphase Flow Handbook. Boca Raton, London, New York: CRC Press Taylor & Francis Group; 2006.
20. Hariprasad DS, Sul B, Liu C, Kiger KT, Altes T, Ruppert K, Reifman J, Wallqvist A. Obstructions in the lower airways lead to altered airflow patterns in the central airway. Respir Physiol Neurobiol. 2020; 272:103311. DOI | PubMed
21. Srivastav VK, Paul AR, Jain A. Capturing the wall turbulence in CFD simulation of human respiratory tract. Mathematics and Computers in Simulation. 2019; 160:23-38. DOI
22. Kim M, Lim J, Kim S, Jee S, Park D. Assessment of the wall-adapting local eddy-viscosity model in transitional boundary layer. Computer Methods in Applied Mechanics and Engineering. 2020; 371:113287. DOI
23. VanDine A, Pham HT, Sarkar S. Investigation of LES models for a stratified shear layer. Computers & Fluids. 2020; 198:104405. DOI
24. Chen J, Gutmark E. Numerical investigation of airflow in an idealized human extra-thoracic airway: a comparison study. Biomech Model Mechanobiol. 2014; 13(1):205-14. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
25. Nicoud F, Ducros F. Subgrid-scale stress modelling based on the square of the velocity gradient tensor. Flow, Turbulence and Combustion. 1999; 62(3):183-200.
26. Lu MZ, Liu Y, Ye JY, Luo HY. Large Eddy simulation of flow in realistic human upper airways with obstructive sleep. Procedia Computer Science. 2014; 29:557-64. DOI
27. Longest PW, Hindle M, Choudhuri SD, Xi J. Comparison of ambient and spray aerosol deposition in a standard induction port and more realistic mouth–throat geometry. Journal of Aerosol Science. 2008; 39(7):572-91. DOI
28. Kiasadegh M, Emdad H, Ahmadi G, Abouali O. Transient numerical simulation of airflow and fibrous particles in a human upper airway model. Journal of Aerosol Science. 2020; 140:105480. DOI
29. Yang XL, Liu Y, So RM, Yang JM. The effect of inlet velocity profile on the bifurcation COPD airway flow. Comput Biol Med. 2006; 36(2):181-94. DOI | PubMed
30. Alzahrany M, Banerjee A. Aerosolized drug delivery in patient-specific lung model during invasive high frequency oscillatory ventilation. Journal of Aerosol Science. 2015; 81:1-20. DOI
31. Comerford A, Förster C, Wall WA. Structured tree impedance outflow boundary conditions for 3D lung simulations. J Biomech Eng. 2010; 132(8):081002. Publisher Full Text | DOI | PubMed [ PMC Free Article ]
32. Jayaraju ST, Brouns M, Lacor C, Belkassem B, Verbanck S. Large eddy and detached eddy simulations of fluid flow and particle deposition in a human mouth–throat. Journal of Aerosol Science. 2008; 39(10):862-75. DOI
33. Chapelier JB, Lodato G, Jameson A. A study on the numerical dissipation of the spectral difference method for freely decaying and wall-bounded turbulence. Computers & Fluids. 2016; 139:261-80. DOI
34. Marati N, Casciola CM, Piva R. Energy cascade and spatial fluxes in wall turbulence. Journal of Fluid Mechanics. 2004; 521:191-215. DOI
35. Aubry N, Guyonnet R, Lima R. Turbulence spectra. Journal of Statistical Physics. 1992; 67(1):203-28. DOI
36. Celik IB, Cehreli ZN, Yavuz I. Index of resolution quality for large eddy simulations. J Fluids Eng. 2005; 127(5):949-58. DOI
37. Kelly JT, Asgharian B, Kimbell JS, Wong BA. Particle deposition in human nasal airway replicas manufactured by different methods. Part I: Inertial regime particles. Aerosol Science and Technology. 2004; 38(11):1063-71.
38. Koullapis PG, Nicolaou L, Kassinos SC. In silico assessment of mouth-throat effects on regional deposition in the upper tracheobronchial airways. Journal of Aerosol Science. 2018; 117:164-88. DOI
39. Krause F, Wenk A, Lacor C, Kreyling WG, Möller W, Verbanck S. Numerical and experimental study on the deposition of nanoparticles in an extrathoracic oral airway model. Journal of Aerosol science. 2013; 57:131-43. DOI
40. Kadota K, Imanaka A, Shimazaki M, Takemiya T, Kubo K, Uchiyama H, Tozuka Y. Effects of inhalation procedure on particle behavior and deposition in the airways analyzed by numerical simulation. Journal of the Taiwan Institute of Chemical Engineers. 2018; 90:44-50. DOI