Prerpints.org logo
Preprint
Article

Influence of Morphological Parameters on the Flow Development Within Human Airways

Altmetrics

Downloads

238

Views

53

Comments

1

A peer-reviewed article of this preprint also exists.

Submitted:

03 February 2023

Posted:

06 February 2023

You are already at the latest version

Alerts
Abstract
Anatomical airways parameters as length, diameter and angles have a strong effect on the flow dynamics. Aiming to explore the effect of variations of the bifurcation angle (BA) and carina rounding radius (CRR) of lower human airways on the respiratory processes, numerical simulations of airflow during inhalation and exhalation were performed using a synthetic bifurcation models. Geometries for the airways models were parameterized based on a set of different BA's and several CRR's. A range of Reynolds numbers (\Reyn) relevant to the human breathing process were selected to analyzed the airflow behaviour. The numerical results show a significant influence of BA and the CRR on the development of the airflow within the airways, and therefore affecting some relevant features of the flow, namely the deformation of velocity profiles, alterations of pressure drop, the flow patterns, and finally enhance or attenuation of wall shear stresses (WSS) appearing during the regular respiratory process. Numerical results show that increases in the bifurcation angle value are accompanied by pressure increases of about 20\%, especially in the regions close to the bifurcation. Similarly, increases in the BA value lead to a reduction in peak shear stresses of up to 70\%. For the ranges of angles and radii explored, an increase in pressure of about 20\% and a reduction in wall shear stress of more than 400\% were obtained by increasing the carina rounding radius. Analysis of the coherent structures and secondary flow patterns also revealed a direct relationship between the location of the vortical structures, the local maxima of the velocity profiles and the local vorticity minima. This relationship was observed for all branches analysed, for both the inhalation and exhalation processes of the respiratory cycle.
Keywords: 
Subject: Engineering  -   Bioengineering

1. Introduction

The human airways, as described by West [52], are divided in three zones: conducting, transitional, and respiratory zones. The conducting zone is intended to transport the air from upper airways to the gas exchange zone, and is composed by the trachea, main bronchus, lobar bronchus and bronchioles. As explained by Minnich & Mathisen [34], the trachea stretches from the lower edge of the cricoid cartilage to the point of bifurcation of the main bronchus, also known as the carina. Some geometrical parameters such as diameters, lengths and angles, among others, make it possible to concretely characterise the types of morphological branches and therefore to establish a hierarchy of generations, as indicated by Singhal et al [43]. The bifurcation angle (BA) is a parameter closely related to lung morphology. Its numerical value depends on a number of factors including chest width [19], and gender and age [2]. Similarly, the choice of techniques and methodologies used to measure the value of the BA has a significant influence on the determination of the BA, with radiography, computed tomography and in-vitro measurements on human cadavers being among the most commonly used [5]. The values of the bifurcation angle reported in the literature vary widely. A range of subcarinal angle variation between 55 ∘ and 65 ∘ for all ages is described by Haskin & Goodman [19], meanwhile in the study carried out by Bipinchandra et al [5] a range between 50 ∘ to 130 ∘ was found. Furthermore, Christou et al [9] reports ranges from 65 . 04 ∘ to 122 . 01 ∘ for men and from 69 . 46 ∘ to 113 . 94 ∘ for women. On the other hand, the study by Sahni et al [42] showed that this parameter has values ranging from 42 ∘ to 75 ∘ in human fetuses. The widening of the bifurcation angle may also be influenced by some physical pathologies such as pericardial fluid accumulation [7], or enlargement of the left atrium [30,36]. These abnormal conditions cause displacement of the bronchi and thus alteration of the bifurcation angle.
Another relevant morphological parameter is the carina rounding radius (CRR). As pointed out by Horsfield et al [23], the rounding radius varies in shape from an almost straight and sharp junction in some cases to a very blunt bifurcation in others. This radius can be related to the diameter of the daughter branch, from which a dimensionless parameter characterising the cross-section as a function of diameter can be estimated. This ratio, although variable, is approximately 0.1 [23]. From a geometrical perspective, Martonen et al [32] describes four potential configurations for this parameter: blunt, parabolic, saddle, and asymmetric. These configurations affect the airflow patterns and conditions within the human airways. For example, asymmetric and saddle configurations are generally linked to the development of larger regions of instability compared to parabolic and symmetric shapes.
The ventilation process and other related phenomena occurring in the human airways during inhalation and exhalation processes have been extensively studied, both numerically and experimentally. The recent comparative study carried out by Xu et al [56] showed how computational fluid dynamics (CFD) techniques adequately capture upper airway fluid phenomena by comparison with experimental results. In fact the research work done by Faizal et al [14] and Tran et al [45] illustrates how these numerical approaches and methods are so robust that they even allow simulating pre and post surgical processes in human airways. The applicability also extends to the study of pathologies such as asthma [8,59], chronic obstructive pulmonary disease -COPD- [18,25], tumors [63], stenosis [35], as well as particle transport phenomena [26,40], coughing and sneezing [20,41], among others. Approaches such as the one developed by Mason et al [33] show through CFD simulations how the aerodynamic stresses on the airways of infants are greater than in adulthood, allowing to characterize some pathologies depending on the age of the patients. Numerical simulation studies generally require the development of a computational domain from a geometrical configuration of the airways. The most commonly used models are the Weibel model [51] and the Horsfield model [23]. For the characterisation of these models it is necessary to define several morphological parameters, among which are the bifurcation angle and the carina rounding radius mentioned above. As described in several works (see [21,24,26,28,46]), the areas near the bifurcation points are of special attention, as they are, for example, associated with particle accumulation and dispersion. In fact, numerical work performed by Balashazy et al [3] to analyze the effect of carina shape on particle behaviour shows how the deposition patterns are quite similar for carina geometries with a straight joint and one with a slight smoothing . Therefore, it is important to characterise the behaviour of fluid flow through the airways, as well as the affectation of flow due to variation of morphological parameters associated with bifurcation, such as BA and CRR. Although a lot of numerical investigations has focused on the analysis of flow through branched systems, the specific effect of BA and CRR on velocity profiles, secondary flow patterns, pressure drops and shear stresses has not been extensively studied.
In order to shed additional light on the effects of morphological parameters on the respiratory cycle, the main objective of the present work was to analyse the effect of variations of the bifurcation angle (BA) and carina rounding radius (CRR) on the respiratory process. The exploration of the effect of BA and CRR was performed by numerical experiments on synthetic models of the human lower airways. Geometric models were used for seven different BA values and three different CRRs. In addition, airflow at different Reynolds numbers was used to explore the influence of these morphological parameters on different flow regimes in the inspiratory and expiratory phases. The comparative analysis performed in this study further deepens the description and characterisation of airflow through the human airways, thus complementing the existing debate on the effects of variations of morphological parameters of the bronchial tree. The results presented in this work focus especially on the description of some flow characteristics that have not been described in depth before, such as wall shear stresses and variability of flow patterns due to morphological changes. This numerical analysis made it possible to explore a large number of configurations that would otherwise involve high experimental complexity.
The present manuscript is organised as follows. The first section of this paper specifies the methodology and the numerical model used in this work, while the second section presents and discusses the results obtained from the numerical simulations for the inhalation (BA and CRR) and exhalation (BA only) processes. Finally, the last section of this manuscript presents the main conclusions and recommendations arising from this work.

2. Methodology

In this study, the effects of variations of two morphological parameters on the breathing process were explored. Bifurcation angle (BA) and carina rounding radius (CRR) were evaluated using independent geometric models and specific boundary conditions. Both morphological parameters were analysed using computational models implemented in OpenFOAM following the assumptions of incompressible and isothermal air as the working fluid. The details of each of these models and their respective methodological aspects are presented below.

2.1. Bifurcation angle (BA). Morphological definition

The analysis of the bifurcation angle was carried out using a Weibel 3-D model of airways, which included the first three generations of human morphology (see [51]). Those generations represent the trachea, main bronchus and secondary bronchus. To simplify the analysis, the notation for airways provided by Christou et al [9] was adopted. The values selected for the geometric parameters of the model are shown in Table 1 and illustrated in Figure 1.
To avoid problems related to flow development in the regions close to the outlets, the length of each of the last generation branches was extended to ensure a minimum length of five times the related outlet diameter. The bifurcation points were considered as straight joints, i.e. the carina rounding radius for these geometries was set equal to zero.
To evaluate the effect of the bifurcation angle (BA), two values of 55 . 9 ∘ and 63 . 1 ∘ were selected from experimental measurements of physiologically relevant subcarinal angles reported by Haskin & Goodman [19]. As shown in Table 2, these angles correspond to the average measured for people between 50 and 60 years old, and the average for people 61 years old or older, respectively. In order to extend the analysis, a model with an additional value for BA of 70 ∘ was also explored, according to the mean angle with minimum resistance discussed by Horsfield & Cumming [22]. Each selected angle value was prescribed as the bifurcation angle for all model generations in each of the numerical experiments.

2.2. Carina rounding radius (CRR). Morphological definition

For the study of the CRR, a synthetic 3-D airways model based on the bifurcation geometry proposed by Lee et al [29] was selected. Parameters and values for the implementation of this model are presented in detail in Table 3 and Figure 2. These generations represent the trachea and main bronchus. As in the BA study, the geometric model was adjusted to ensure a minimum length of the branches of the last generation of at least five times the related outflow diameter. Specifically, the parameters R d and α describe the curvature of the bifurcation. The CRR was non-dimensionalised using the following expression,
r c = C R R r t
where r c is the dimensionless rounding radius and r t the trachea radius. As reported by Horsfield et al [23] and Balashazy et al [3], for human airways the average value of r c ranges around 0.1 . In fact, dimensionless rounding radius of 0.07 and 0.14 for symmetric and asymmetric airways, respectively, were used in the lung architecture characterisation work carried out by Lee et al [29]. Based on the above, three dimensionless carina radii, equivalent to 0, 0.07 and 0.14 , were selected in the present work.

2.3. Governing equations and computational model.

The ventilation process within human beings takes place at relatively low Reynolds numbers regimes, and at very low Mach numbers, even for agitated breathing situations ([31,32]). Therefore, for both BA and CRR experiments, the fluid flow through the synthetic human airways was assumed to be laminar and incompressible. Additionally, isothermal air was assumed as the working fluid, in line with the predefined flow conditions. The governing equations selected for this fluid flow regime were the conservation of mass equation and the Navier-Stokes (conservation of momentum) equations, which, for incompressible and isothermal flow, can be formulated as,
∇ · v → = 0
∂ u i ∂ t + ∇ · ( u i v → ) = − 1 ρ ∇ p + ∇ · ( ν ∇ u i )
where u i and v → are the i-th velocity component and the velocity vector field, respectively Equally, in Equation (3), p is the fluid pressure and ν = μ / ρ its kinematic viscosity with ρ standing for the density of the fluid. For the purposes of this work, the Reynolds number has been defined as
Re = U a v e D ν
Where U a v e = 0.5 U m a x . Air properties were defined at ambient temperature of 15o, with ρ = 1.23 [ k g / m 3 ] and ν = 1.48 e − 5 [ m 2 / s ] .
The two breathing processes, i.e., inhalation and exhalation, were decoupled and studied as independent processes. The main objective of this work was to characterize and establish general patterns of airflow within a portion of the human airways. By simulating the processes independently we were able to guarantee the stabilization of the flow and thus the stabilization of the patterns reported in the present manuscript, mainly aiming at getting a general description and overview of the flow. It is also noteworthy that a large body of research in this field has been performed employing this consideration of independent processes, including experimental and numerical studies (see [10,15,57,58,60,61,62]). For instance, Zhao & Lieber [61,62] characterized the flow through the airways experimentally for the inhalation and exhalation states independently; Yang et al [57] characterized airway inflow patterns as a single process, and Zhang et al [60] described the effect of tracheal curvature on respiratory inhalation as an independent stage. A large number of numerical studies on particle dispersion, such as those developed by Yu et al [58] and Comer et al [10], consider a single breathing state to characterize the phenomena and obtain patterns of behavior. Accordingly, considering the large portion of published work on the modeling of the human respiratory process based on the one-way respiration flow assumption, instead of a more realistic oscillatory flow, and the success of this approach in modeling a number of different airflow conditions, the same methodology was used to explore the cases considered in the present work. To ensure proper flow development in the inhalation process simulations, a parabolic velocity profile was adopted. This profile was prescribed as symmetrical around the inlet axis, i.e. the trachea, and is given by:
U ( r ) = U m a x 1 − r R 2
Where U ( r ) is the inlet velocity, U m a x is the maximun velocity, r a radial position and R the inlet radius. Human breathing rates range from 0.2 L / s to 2.5 L / s [31,32], which corresponds to a range of Reynolds numbers (Re) between 200 and 2800. This range was observed by Adler & Brücker [1] in an experimental investigation using a pulsatile flow, where a maximum value of Re = 3000 was reached. Different Reynolds number values were explored in the present study, however, for brevity only the results for two values, representative of the extremes of the breathing process and corresponding to Re = 500 and Re = 2000 , are presented. At the Outlets, a uniform pressure boundary condition was prescribed for all outlets. As indicated by Bauer & Brücker [4], this consideration allows the flow to adjust to different pressure gradients across different airway generations. Therefore, the velocity condition at the outlets is given by
∂ v → ∂ n = 0
In the surface of the wall a Dirichlet no-slip boundary condition was applied to the velocity field at all airway walls,
v → = 0
The exhalation study was performed only for BA analysis. While in inhalation stage the inlet was the trachea and the outlets were the bronchi, for exhalation stage the bronchi were the inlets and the trachea was the outlet. A parabolic velocity profiles (see Equation (5)) were prescribed at each of the inlets. These profiles were defined symmetric with respect to the axial lines of each last generation branch, and their maximum value is equal to the maximum velocity obtained at the outlets in the respective inhalation cases. At the Outlet a uniform pressure condition was set up. The velocity condition is described by Equation (6). Newly a Dirichlet no-slip boundary condition was applied to the velocity field at all airway walls, as described in Equation (7).
All boundaries were considered rigid. Investigations such as those carried out by Wall & Rabczuk [47], Wang et al [48] and Koombua et al [27] show that under normal breathing conditions, i.e., low Reynolds numbers and healthy airway geometries, the wall deformation is relatively small, and therefore the effect on airflow is not highly relevant. The study performed by Wall & Rabczuk [47] shows that the difference between the numerical results obtained with fluid-structure interaction simulations versus CFD results was minimal, and in fact the flow fields were very similar for small Tidal breathing volumes, i.e., for considerations close to normal breathing in physiological conditions.
The governing equations were solved numerically using the open source computational fluid dynamics (CFD) suite OpenFOAM, which offers different schemes and numerical methods to solve different fluid flow phenomena. For the selected flow conditions and mathematical model, a finite volume method (FVM) approach was adopted. OpenFOAM (OF) offers a number of the so-called solvers and for the purposes of the present study the pisoFoam solver was selected, which is suitable for unsteady and incompressible laminar flows. The simulation setup was configured to maintain a Courant number less than 1, with a time step size equal to 1 × 10 − 4 . This consideration ensures that the solution has a stable behavior during the simulation time [13,49,53]. Temporal integration was performed by using a backward scheme and the spatial discretisation was solved using second-order interpolation schemes for the gradient and divergence terms. The Geometric Agglomerated Algebraic Multigrid Solver (GAMG) was used to solve the Poisson-pressure equations, whereas the velocity field was obtained using the smoothSolver scheme available in OF. Pressure-velocity coupling was accomplished by using the PISO algorithm. All the numerical simulations of this study were performed in parallel using a Simple Domain Decomposition. The unstructured mesh generator snappyHexMesh was employed to generate the computational grid for the different geometries covered in this work. An example of one of the finite volumen mesh for BA’s study is illustrated in Figure 3a. As can be seen in the cross section of the branch, a finer refinement was implemented near the walls. The average number of tetrahedral cells in the different computational models was 3.5 × 10 6 . The finite volumen meshes for CRR’s study are shown in Figure 3b. In order to obtain a zone of finer cells near the carina, particular sub-regions of refinement were defined. For this cases, the average number of tetrahedral cells was 1.5 × 10 6 .
To ensure that the resolution of the mesh was capable of capturing the associated flow phenomena, the computational mesh was configured so that the wall-normal distance of the cell centres neighbouring the airway walls was less than the resolution threshold of the viscous sublayer. This value was expressed in non-dimensional form in terms of wall-units, generically referred to as y + . This parameter can be evaluated as,
y + = u τ y ν
where ν is the fluid’s kinematic viscosity, u τ is the friction velocity defined as
u τ = τ w ρ
the variable y is the wall-normal distance from the cell centre to the bounding wall, and τ w the local wall shear stress. As the larger velocities appear in the trachea, y + was computed in that branch, as shown in Figure 4. The parameter τ w was calculated at the maximun velocity of our test set, i.e. Re = 2000 . The average value obtained by applying Equations (8) and (9) was y + = 0.258 . Considering that this value is less than 1, it is possible to conclude that the meshes constructed and employed in the present work were fine enough to accurately capture the flow phenomenon for the different regimes considered here.

2.4. Validation of the numerical model.

The numerical model applied in this study was tested and validated against experimental results obtained by Zhao & Lieber for both inhalation [61] and exhalation [62] flow. The model selected for this validation is the simple bifurcation geometry developed for the CRR study (see Figure 2), since this configuration is similar to the geometry of the experimental study, with an opening angle of 70° and with a rounding radius of 0. The comparison was performed using the normalized profiles with the maximum velocity U m a x for Re = 500 . The comparative results for flow in inhalation are illustrated in Figure 5a and Figure 5b. These profiles are taken at the end of the second generation branch. It can be seen how the main flow features are accurately captured, particularly the narrowing and deformation of the velocity profile towards the inner wall of the branch (on the coronal plane), and the formation of the M-shape of the profile (on the sagittal plane). A similar result is obtained in the comparison of the results in the exhalation stage. As evidenced in Figure 5c and Figure 5d, both the peak velocity centered on the axial branch axis in the coronal plane and the M-shaped profile over the sagittal plane are faithfully captured. These comparative results confirm that the numerical model selected for the exploration of the effect of BA and CRR on airflow was accurate and suitable for the purposes of our study.
Our discussion and results for both the effects of BA and CRR variations focus on the analysis of velocity profiles, secondary flow patterns, pressure drops and wall shear stresses (WSS). These indicators are generally analysed in the coronal and sagittal planes. In the present work, these values were explored specifically in the areas and at the points shown in Figure 6.

3. Results I. Effect of variations of the bifurcation angle (BA)

3.1. Velocity Profiles

One of the main features of any flow pattern is the mean velocity profile. For the intended analysis, mean velocity profiles on the coronal plane for the selected BAs were constructed and are shown in Figure 7 for flows at Re = 500 and Re = 2000 . These curves were built using a normalized length d / D which was measured at the centre of each branch: Left Main Bronchus (LMB), Left Upper Lobar Bronchus (LUB) and Left Lower Lobar Bronchus (LLB), following the nomenclature described in Figure 1. The unit of the magnitude of the mean velocity is [m/s]. Velocity averages were performed only when the flow was fully developed, i.e. in the steady state of the simulation. As reported by different authors (see [6,31,44,54,61]), a narrowing of the velocity profile, as well as a shift of the maximum peak from the outer to the inner walls, are among the main effects of flow through a bifurcation. Observing Figure 7, it can be seen how these alterations of the velocity profiles are affected by changes in the bifurcation angle, especially at high Re numbers. For the Re = 500 regime the velocity profiles measured at the three selected positions (LMB, LUB and LLB) present a similar behavior, with a single velocity peak tending towards the inner wall of the branch, as shown in Figure 7a to Figure 7c. Noteworthy, at the lower Re regime, the BA mildly affects the narrowing of the profile, where there is a slightly increase of the peak value at lower BA values, although this effect is almost imperceptible in the LUB branch. In turn, for Re = 2000 , only in the LMB branch, Figure 7d, the behavior described above is preserved. For the LUB branch, shown in Figure 7e, a two-peaks profile is obtained for the larger angles. In the case of the LLB branch, Figure 7f, a deceleration in the center of the branch leads to the development of an even clearer two-peaks profile, an effect that becomes more noticeable as the BA decreases. It seems clear then that the BA deformation effect on the velocity profiles, at least in comparison with symmetrical profile patterns, is stronger as the Reynolds number increases.
Figure 8 show the velocity profiles taken in the sagittal plane for the selected angles at Re = 500 and Re = 2000 . These profiles show mirror symmetry around the axis of the branches, and a characteristic "M" shape can be identified. This shape was described by Schroter & Sudlow [44] and it illustrates how the flow velocity increases near the walls and decreases towards the branch axis. As shown in Figure 8a and Figure 8c for Re = 500 , and Figure 8d and Figure 8f for Re = 2000 , this effect is more noticeable for the LMB and LLB branches as BA increases. In contrast, for the LUB branch, Figure 8b and Figure 8e, it is more pronounced as BA decreases. Again, it is evident that the velocity profiles taken in the sagittal plane are more affected by the BA as the Reynolds number becomes larger.
Results for the exhalation stage, presented in Figure 9, show how the velocity profiles are symmetrical both in coronal and sagittal planes, a result that was previously observed and reported by Corieri [11], Schroter & Sudlow [44]. The existence of a single velocity peak centered in the middle of the branch in the coronal plane is illustrated in Figure 9a. In turn, the velocity profiles in the sagittal plane develop again the M-shape, as shown in Figure 9b, which had been obtained for the branches in the inhalation stage. These effects are observed for Re = 500 and are more noticeable for lower BA values. The velocity profiles in the coronal and sagittal planes for Re = 2000 are shown in Figure 9c and Figure 9d, respectively. These plots show how a higher flow velocity leads to the generation of three peaks located both in the center and near the branch walls. For the profiles taken on the coronal plane this effect is more pronounced for larger BAs, contrary to the sagittal plane where the effect increases with lower BA’s.

3.2. Flow patterns

The effect of the changes of the selected morphological parameters was also explored by studying coherent vortical structures. In particular, Figure 10 shows the behavior of the vortical structures through the branched structure by applying the Q-criterion, for both Reynolds numbers at the inhalation stage. This criterion is based on the second invariant of the velocity gradient tensor and computed as
Q = 1 2 ( | | Ω | | 2 − | | S | | 2 )
where Ω is the rotation tensor and S is the rate of strain tensor. The dimensionless factor Q * = Q ( D 2 ) / ( U b u l k ) 2 described by Han et al [17] was adopted. Likewise, the color scales for the structures are presented using the dimensionless velocity U * = U / U b u l k . The two vortex cores generated after each bifurcation point for all branches are shown in Figure 10a for Re = 500 . These structures are symmetric on the coronal plane and their location and quantity are changing with respect to the number of Reynolds, as evidenced by the isocontours for Re = 2000 in Figure 10b. The presence of four vortices for the last generation (LLB and LUB branches) is evidenced. However, it is striking that the patterns obtained using the Q-criterion show no significant changes for the range of bifurcation angles selected in this study. As illustrated in Figure 11, the vortical structures for an angle of B A = 55 . 9 ∘ are similar to those for B A = 70 ∘ at Re = 2000 .
A deeper analysis of these structures was carried out through the study of secondary flow patterns (SFP’s). This analysis was performed using a line integral convolution projection, commonly known as Surface LIC. Using the surface LIC method it is possible to establish visual analogies with temporal stability points, and to identify spatial points that act as orbit centres or as saddle points. In the present work, a correlation between these SFPs and velocity and vorticity profiles was obtained for all the bifurcation angles studied. By performing a juxtaposition of the Surface LIC plots with the velocity and vorticity profiles, a particular relationship between the vortical structures and the profiles was observed. An example of this characteristic link for the inhalation stage is presented in Figure 12 for a model featuring a BA equals to 55 . 9 ∘ . Examining the secondary flow patterns it seems clear that in all branches dean vortices [12] and saddle points are developed. These vortical structures are symmetric with respect to the coronal plane and their quantity and location depend on the Reynolds number, as discussed above with the coherent structures. If the velocity and vorticity profiles are superposed to each other using a non-dimensional axis, a composed profile is obtained as those presented in Figure 12. In the present work such profiles are called “V-W" profiles. By examining these profiles, it is possible to observe how the velocity peak coincides with the point of minimum vorticity, and which also consistently coincides with the location of the saddle point on the secondary flow. As illustrated in Figure 12a and Figure 12b, the V-W curves for Re = 500 present a single velocity peak and a vorticity minimum which are aligned with the saddle point. For higher flow velocities two new vortices are generated and the saddle point is relocated. As previously shown in Figure 7f, for Re = 2000 there are two velocity peaks in the LLB branch. By analyzing Figure 12c, it is observed how the peak of higher magnitude is also aligned with a local minimum of vorticity and with the central saddle point. In turn, the secondary velocity peak is also aligned with another local minimum of vorticity and with another not so well defined saddle point.
The coherent structures for the exhalation stage at different Reynolds numbers are shown in Figure 13. Similarly to the structures observed in the inhalation case, the vortex cores after the first bifurcation point give an indication of the existence of four vortices in the LMB branch. In the trachea, from four to eight vortices are developed as the Reynolds number increases up to Re = 2000 . The comparison of these structures for B A = 55 . 9 ∘ and B A = 70 ∘ is illustrated in Figure 14. A slight difference between the vortex core lengths in the LBM branch is observed. A longer vortex core length is evident for the larger bifurcation angle.
The secondary flow patters and the V-W profiles at exhalation stage are shown in Figure 15. These flow patterns and vorticity profiles are symmetric with respect to the coronal and sagittal planes. Figure 15a evidences the existance of four vortices rotating around the axial axis of the trachea for Re = 500 , in line with the description derived from coherent structures. The saddle point is located exactly at the center of the branch, as well as the peak velocity and the local minimum vorticity. Similarly, for the LMB branch at Re = 500 the saddle point remains centered on the axial axis and aligned with the peak velocity and minimum vorticity, as illustrates Figure 15c. As described in Figure 9d, the velocity profiles over the sagittal plane for Re = 2000 show three peaks. Figure 15b illustrates how these peaks are coincident with three local minima of vorticity and three saddle points visible through the flow patterns. For this case, eight vortices were clearly identifiable.

3.3. Pressure behaviour

Pressure drop behaviour, as explained by Pedley et al [38], is directly influenced by changes in kinetic energy and by the dissipation of viscous energy. Geometry configuration and the definition of parameters such as length and cross section impact the pressure behavior in branched systems ([39]). In fact, the numerical study in micro-channel networks developed by Wang et al [50] illustrates how the growth of the bifurcation angle influences the increase of pressure drops for branched systems of rectangular cross section. Plots of Δ P versus a normalizad length along the branches axial axes are presented in Figure 16. The unit of the magnitude of the pressure is [Pa]. In all generations for all the cases explored a local pressure growth at the beginning of the axial line is evidenced. This behavior is a result of the increase of pressure gradients generated by the bifurcation points. As might be expected, a higher magnitude in pressure drop is the result of an increase in BA for all branches. This local growth is slightly noticeable for the LLB branch (Figure 16c and Figure 16f) with respect to the LUB branch (Figure 16b and Figure 16e). This difference is related to the velocity profiles over the coronal and sagittal plane shown in Figure 7 and Figure 8, respectively. For the LLB branch, the velocity magnitudes are slightly larger than for the LUB branch. The maximum of this local pressure increase between the angle of B A = 55 . 9 ∘ and B A = 70 ∘ is up to 20% for the first bifurcation point.
The pressure drop profiles in the exhalation state for Re = 500 are shown in Figure 17. these profiles are plotted in the upstream direction, i.e., from the last generation towards the trachea. Contrary to the results obtained in the inhalation state, no particular behavior derived from the bifurcation angles is identified for the exhalation stage.

3.4. Wall Shear Stresses

During inhalation and exhalation stages, the maximum shear stresses developed are located on the inner and externals wall, respectively, as reported by different authors [16,37,55]. In a similar manner as discussed with the pressure drop, and as it might be expected, this behaviour of the WSS is also strongly related to the specific changes of the velocity profiles reported in this paper. Particularly, the previously discussed tendency of displacement of the peak of the profiles (shown earlier in Figure 7) is always associated to the appearance of regions of maximum wall shear stress.
Wall shear stress was measured along to the inner wall (for inhalation) and in a perpendicular direction of the inner wall (for exhalation) on each of the branches. In the different plots (Figure 18 and Figure 19), the length is presented in normalized form ( x / L , where L is the total length of the branch). In all cases, for the x-label, the value of "0" and "1" refers to the upstream and downstream, respectively. The unit of the magnitude of the WSS is [Pa]. By analysing the effect of the bifurcation angle (BA) on the wall shear stress distribution for the inhalation process, it is possible to determine that, as observed in different studies, the maximum shear stress is located towards the intersection of the bifurcation, that is in the region nearby the joints. It is further observed that for the different branches, as the BA considered is augmented, the magnitude of the maximum WSS decreases, giving as a result that, for instance, the highest figures for WSS are observed for the smallest angle, which in the present study is B A = 55 . 9 ∘ . The decrease of the WSS between the angle of B A = 55 . 9 ∘ and B A = 70 ∘ is up to 70% at the the bifurcation point. Examining the WSS in each of the branches, along the local axial direction and therefore progressing through the branching line, the WSS seem to collapse reaching an average common value regardless of the BA explored. This trend, however, is sightly disrupted in LLB (Figure 18c and Figure 18f), where it is possible to appreciate a very subtle difference between the Mean values of wall shear stress attained for B A = 63 . 1 ∘ in comparison to the other BA explored. Thus, it can be established that the shear forces vary with respect to the bifurcation angle only in the zones or regions near the bifurcation point, but converging to an average value downstream of such a point when considering inhalation.
The results obtained for the exhalation process, and presented in Figure 19, show a complete different picture for the distribution of the wall shear stresses, at least in comparison with the inhalation process described just above. For instance, it is possible to observe that exhalation brings about an increase of the WSS in the first half of each branch, so the maximum values of WSS are present between 0.2 < x / L < 0.5 . In this case x / L is measured from the upstream bifurcation following the downstream direction and, as mentioned previously, using as normalising factor the length of the respective branch L. In any case, for all the branches, there is a reduction of the WSS that seems to have a similar behaviour for the different angles considered in this work. In this exhalation process, however, a converging or collapsing trend to an average value is not clear. Noteworthy, the maximum WSS in each branch is accompanied by strong oscillations, albeit this oscillatory tendency is stronger for LMB than for the TRA. This particular effect may be related to the onset of flow instabilities and eventually some perturbations which grow in the flow direction as a result of the increment of the local flow rate and cross section area. A clear indication of this phenomena is increase of the local Reynolds number as the fluid flows downstream towards the TRA.

4. Results II. Effect of variations of the carina rounding radius (CRR)

4.1. Velocity profiles

Mean velocity profiles on the coronal plane for the selected CRRs were constructed and are shown in Figure 20 and Figure 21 for flows at Re = 500 and Re = 2000 . Newly, these curves were built using a normalized length d / D and these were measured in sections at 20% and 50% of the total LMB branch length, with the main intention of analyzing the velocity profiles both at a point close to the carina and in the middle of the branch. Figure 20 shows how the rounding radius has a slight effect on the velocity profiles measured at 20% of the LMB length. The narrowing of the profile toward the inner face of the branch is more noticeable as the CRR is larger. As the flow is transported downstream this behavior becomes almost imperceptible, as can be seen in Figure 21 for velocity profiles measured at 50% of the LMB length.
Figure 22 and Figure 23 show the velocity profiles over the sagittal plane for Re = 500 and Re = 2000 . Theses measurement were taken at 20% and 50% of the LMB branch length. In contrast to the coronal plane profiles, the effect of the CRR variation is clearly noticeable. The flow experiences a greater acceleration as the CRR is larger for both Re = 500 and Re = 2000 , respectively. It is also possible to observe how the profiles near the bifurcation develop the maximum peak at the center of the branch, as illustrated in Figure 22. The profiles take on the characteristic M-shape previously described in the results for the BA as the flow moves downstream, as shown in Figure 23.

4.2. Flow patterns

The coherent structures for the different Reynolds numbers are shown in Figure 24. These structures are calculated through the Q-criterion, with the same procedure described in the BA section. These isocontours are provided with the color scale of the dimensionless velocity U * . For this simple bifurcation geometry it is possible to appreciate the existence of two vortex cores after the bifurcation, with a symmetric behaviour with respect to the coronal plane. Figure 25 illustrates the comparison of coherent structures for two rounding radius explored in this investigation. As can be seen, no effect derived from the variation of the CRR is clearly visible since the structures are very similar to each other.
Vorticity profiles on the coronal plane for different for the selected CRR are shown in Figure 26 for flows at Re = 500 and Re = 2000 . These curves were built using a normalized length d/D and a normalized vorticity ω / ω m a x and were taken at 20% of the total LMB branch length. The vorticity measured at 50% branch length is not shown in this paper since it does not show any singular behavioral pattern concerning the CRR variation. The maximum vorticity intensity is located on the inner wall of the branch for Re = 500 , as shown in Figure 26a. The existence of a local vorticity peak located near the center of the branch is also evidenced. The increase of the magnitude of this local peak is associated with the increase of the CRR. This effect is preserved as the flow accelerates, as seen in Figure 26b for Re = 2000 . This behavior indicates that the rounding radius influences to some extent the rotation that the flow experiences as it passes through the bifurcation.
As well as in the BA study, the secondary flow patterns were performed using the Surface LIC method. Figure 27 shows the juxtaposition of the Surface LIC plots with the velocity and vorticity profiles in the sections at 20% and 50% of the branch length. These plots corresponding to the geometry with r c = 0 at Re = 2000 . Two dean vortices and a saddle point are clearly identified, in line with the coherent structures described above. The relationship of these vortical structures to the V-W curves is again observed. As illustrated in Figure 27a, at 20% of the branch length the vorticity has not reached its global minimum. This value is reached as the flow advances downstream and its location coincides with the peak velocity, as shown in Figure 27b.

4.3. Pressure behaviour

The pressure drop across the LMB branch for the selected RRCs is shown in Figure 28 for flows at Re = 500 and Re = 2000 . These curves were built using a normalized length x/L which was measured over the branch inner wall. The plots show only 20% of the length since beyond this point the curves do not show any CRR-derived variation. The maximum peak pressure is located exactly at the bifurcation point and is higher as the CRR increases. At the same time it can be observed that the pressure drop is also higher for larger CRRs and smaller for a straight junction. This effect is evident in the two Reynolds numbers analyzed. The maximum of this local pressure increase between the radius of r c = 0 and r c = 0.14 is up to 20% for the bifurcation point.

4.4. Wall Shear Stresses

As with the pressure drop, the wall shear stress behavior is measured on the internal wall of the branch and expressed using the normalized length x/L. These curves are plotted up to 50% of the length since after this point the effects are negligible. Figure 29a shows how the maximum WSS for Re = 500 is developed by the bifurcation with straight joint r c = 0 . As expected, this point acts as a kind of stress concentrator. The location of the WSS maximum moves away from the carina in the downstream direction as the CRR increases and the joint becomes smoother. This displacement tendency of the maximum WSS is also observed for Re = 2000 , as illustrated in Figure 29b. The decrease of the WSS between the rounding radius of r c = 0 and r c = 0.14 is up to 400 % at the bifurcation point. Although the WSS is still higher for the straight joint at the exact point of the carina, for this case the maximum WSS is developed by the geometry with radius r c = 0.14 . These results are in line with the velocity profiles in the sagittal plane. The increase of the rounding radius generates a higher flow acceleration and therefore a higher accumulation of wall shear stresses.

5. Discussion and Conclusions

The present numerical study aimed to produce an assessment of the effect of bifurcation angle (BA) and the carina rounding radius (CRR) on the behaviour of the flow through synthetic respiratory airway models. The main focus was to detail and quantify the influence of these morphological parameters on the respiratory process examining both the inhalation and exhalation stages. Numerical experiments were performed on airway models set up for three different angles and three rounding radius, using two symmetrical models (Weibel et al [51] for BA analysis and Lee et al [29] model for CRR). The respiratory mechanisms of human inhalation and exhalation were decoupled and hence considered as independent processes. Analysis of the results derived from the variations of the bifurcation angle indicate a remarkable influence of this morphological parameter on the flow behavior through the airways. For the inhaled state, the narrowing of the velocity profile in the coronal plane with a tendency towards the inner wall of the airway becomes more pronounced as the BA increases. In turn, this variation in BA drives the acceleration near the walls of the branches in the profiles measured in the sagittal plane. A particular relationship between the vortical structures and the velocity an vorticity profiles was identified: the coincidence in the location of the peak velocity profile, the miminum vorticity magnitud and the saddle point seen through the secondary flow patterns. This behavior is observed for the different flow velocities analyzed. Moreover, a local pressure growth is generated by gradients at the bifurcation points. This local growth has a greater magnitude for the larger BAs, with a magnitude difference of approximately 20% between the largest and smallest angles. At the same time, it is noteworthy that this larger angles impacts on the decrease of wall shear stresses developed near the bifurcation, reducing the maximum shear value between explored angles by almost 70%. The variation of the BA also has a noticeable effect on the velocity profiles in the exhalation stage. The mirror symmetry around the branch axis in both the coronal and sagittal planes is the most notable feature in these profiles. As the flow velocity increases from Re = 500 to Re = 2000 an acceleration near the branch walls is evident. This acceleration is greater for the lower BAs in the sagittal plane and for the higher BAs in the coronal plane. Wall shear stress for this stage exhibits a particular increase in the first half of the branch for all cases. This behaviour may be related to the resistance and oscillations produced by the changes in the cross section area as the flow is transported through the branches. The relationship between the V-W profiles and the vortical structures is also observed at exhalation. In fact, this coincidence in the location of the vortical structures is present in all existing velocity peaks. Contrary to the results obtained at inhalation stage, the effect of the BA varations on the behaviour of the pressure flow was almost unnoticeable.
As with the BA, the results obtained in the CRR analysis show a noticeable impact on the flow. The increase in the CRR causes a greater acceleration of the flow observed in the velocity profiles over the sagittal plane. In turn, it also slightly influences the narrowing of the profile toward the inner walls of the airways measured in the coronal plane. The variation of this morphological parameter also has an impact on the flow rotation. As the radius increases, the rotation observed through the vorticity curves on the coronal plane is greater. The relationship of flow patterns to velocity and vorticity profiles is also evidenced by variations in the CRR. The carina is the point of interest for the analysis of pressure drop and wall shear stress. When the bifurcation is considered as a straight junction (with CRR equal to 0) the maximum pressure at the carina decreases and the WSS is maximum. With the increase of the CRR, i.e. the smoothing of the junction, the pressure drop increases approximately 20% between the largest and smallest CRR, and the maximum WSS is reduced by up to 400%. In turn, the maximum shear stress moves away downstream of the carina (exact point of attachment) as the rounding radius becomes larger.

Author Contributions

Conceptualization, AS Espinosa-Moreno and CA Duque-Daza; Formal analysis, AS Espinosa-Moreno and CA Duque-Daza; Investigation, AS Espinosa-Moreno; Methodology, AS Espinosa-Moreno; Resources, CA Duque-Daza and DA Garzón-Alvarado; Software, AS Espinosa-Moreno; Supervision, CA Duque-Daza and DA Garzón-Alvarado; Validation, AS Espinosa-Moreno and CA Duque-Daza; Writing – original draft, AS Espinosa-Moreno; Writing – review & editing, CA Duque-Daza and DA Garzón-Alvarado.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

The work was technical supported by the investigation group GNUM (Grupo de Modelado y Metodos Numericos en Ingenieria) of the Universidad Nacional de Colombia.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature and Abbreviations

Nomenclature

The following nomenclature are used in this manuscript:
D Branch diameter [ m m ]
R Branch radius [ m m ]
L Branch length [ m m ]
r c Dimensionless rounding radius
y Grid length in the normal direction of the wall [m]
U Inlet velocity [ m / s ]]
U m a x Maximun velocity [ m / s ]
U a v e Average velocity [ m / s ]
u τ Friction velocity [ m / s ]
p Fluid pressure [ P a ]
R e Reynolds number
ρ fluid density [ k g / m 3 ]
ν kinematic viscosity [ m 2 / s ]
μ Dynamic viscosity [ P a . s ]

Abbreviations

The following abbreviations are used in this manuscript:
CFD Computational Fluids Dynamics
WSS Wall Shear Stress
BA Bifurcation angle
CRR Carina rounding radius
TRA Trachea
RMB Right Main Bronchus
LMB Left Main Bronchus
RUB Right Upper Lobar Bronchus
LUB Left Upper Lobar Bronchus
RIB Right Intermediate Bronchus
LLB Left Lower Lobar Bronchus

References

  1. Adler, K. and Brücker, C. Dynamic flow in a realistic model of the upper human lung airways. Experiments in Fluids 2007, 43(2-3), 411.
  2. Alavi, S. M., Keats, T. E., & O’Brien, W. M. The angle of tracheal bifurcation: its normal mensuration. American Journal of Roentgenology 1970, 108(3), 546–549.
  3. Balashazy, I., Heistracher, T., & HOFMANN, W. (1996). Air flow and particle deposition patterns in bronchial airway bifurcations: the effect of different CFD models and bifurcation geometries. Journal of Aerosol Medicine 1996, 9(3), 287-301.
  4. Bauer, K. and Brücker, C.The influence of airway tree geometry and ventilation frequency on airflow distribution. Journal of biomechanical engineering 2015, 137(8) 081001.
  5. Bipinchandra, K., Waheed, A. R. A., Yadav, N., & Diwan, C. Study of sub carinal angle of human trachea by computerized tomography. International Journal of Anatomy and Researchy 2016, 4(3), 2828–2832.
  6. Calay, R., Kurujareon, J., & Holdø, A. E. Numerical simulation of respiratory flow patterns within human lung. Respiratory physiology & neurobiology 2002, 130(2), 201–221.
  7. Chen, J., Putman, C. E., Hedlund, L. W., Dahmash, N., & Roberts, L. Widening of the subcarinal angle by pericardial effusion. American Journal of Roentgenology 1982, 139(5) 883–887.
  8. Choi, S., Yoon, S., Jeon, J., Zou, C., Choi, J., Tawhai, M. H., & Lin, C. L. 1D network simulations for evaluating regional flow and pressure distributions in healthy and asthmatic human lungs. Journal of Applied Physiology 2019, 127(1), 122-133.
  9. Christou, S., Chatziathanasiou, T., Angeli, S., Koullapis, P., Stylianou, F., Sznitman, J. & Kassinos, S. C. Anatomical variability in the upper tracheobronchial tree: sex-based differences and implications for personalized inhalation therapies. Journal of Applied Physiology, 2021, 130(3), 678-707.
  10. Comer, J. K., Kleinstreuer, C., & Kim, C. S. Flow structures and particle deposition patterns in double-bifurcation airway models. Part 2. Aerosol transport and deposition. Journal of Fluid Mechanics 2001, 435, 55-80.
  11. Corieri, P. Experimental and numerical investigation of flows in bifurcations within lung airways. Doctoral dissertation, Ph. D. thesis, von Karman Institute for Fluid Dynamics, Université Libre de Bruxelles and Rheinisch-Westfälische Technische Hochschule Aachen, 1994.
  12. Dean, W. R. Note on the motion of fluid in a curved pipe. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 1927, 4(20), 208–223. [Google Scholar] [CrossRef]
  13. Duque-Daza, C. A., Ramirez-Pastran, J., & Lain, S. Influence of particle mass fraction over the turbulent behaviour of an incompressible particle-laden flow. Fluids 2020, 6(11), 374.
  14. Faizal, W. M., Ghazali, N. N. N., Khor, C. Y., Badruddin, I. A., Zainon, M. Z., Yazid, A. A., & Razi, R. M. Computational fluid dynamics modelling of human upper airway: A review. Computer methods and programs in biomedicine 2020, 196, 105627.
  15. Fresconi, F. E. & Prasad, A. K. Secondary velocity fields in the conducting airways of the human lung. Journal of biomechanical engineering 2007, 129(5), 722–732.
  16. Green, A. S. Modelling of peak-flow wall shear stress in major airways of the lung. Journal of Biomechanics, 2004; 37, 661–667. [Google Scholar]
  17. Han, F., Liu, Y., Lan, Q., Li, W., & Wang, Z. CFD investigation on secondary flow characteristics in double-curved subsea pipelines with different spatial structures. Journal of Marine Science and Engineering 2022, 10(9), 1264.
  18. Hariprasad, D. S., Sul, B., Liu, C., Kiger, K. T., Altes, T., Ruppert, K., & Wallqvist, A. Obstructions in the lower airways lead to altered airflow patterns in the central airway. Respiratory Physiology & Neurobiology 2020, 272, 103311.
  19. Haskin, P. H. & Goodman, L. Normal tracheal bifurcation angle: a reassessment. American Journal of Roentgenology 1982, 139(5), 879–882.
  20. Hassani, K., & Khorramymehr, S. In silico investigation of sneezing in a full real human upper airway using computational fluid dynamics method. Computer methods and programs in biomedicine 2019, 177, 203–209.
  21. Hofmann, W. (2011). Modelling inhaled particle deposition in the human lung: a review. Journal of Aerosol Science 2011, 42(10), 693–724.
  22. Horsfield, K. and Cumming, G. Angles of branching and diameters of branches in the human bronchial tree. The Bulletin of mathematical biophysics 1967, 29(2), 245–259. [CrossRef] [PubMed]
  23. Horsfield, K., Dart, G., Olson, D. E., Filley, G. F.and Cumming, G. Models of the human bronchial tree. Journal of applied physiology 1971, 31, 207–217.
  24. Islam, M. S., Paul, G., Ong, H. X., Young, P. M., Gu, Y. T., & Saha, S. C. A review of respiratory anatomical development, air flow characterization and particle deposition. International journal of environmental research and public health 2020, 17(2), 380.
  25. Kadota, K., Matsumoto, K., Uchiyama, H., Tobita, S., Maeda, M., Maki, D., & Tozuka, Y. In silico evaluation of particle transport and deposition in the airways of individual patients with chronic obstructive pulmonary disease. European Journal of Pharmaceutics and Biopharmaceutics 2022, 174, 10–19.
  26. Kim, Y., Tong, Z., Chan, H., & Yang, R. CFD modelling of air and particle flows in different airway models. Journal of Aerosol Science 2019, 134, 14–28.
  27. Koombua, K., Pidaparti, R. M., Longest, P. W., & Ward, K. R. Computational analysis of fluid characteristics in rigid and flexible human respiratory airway models. Engineering Applications of Computational Fluid Mechanics 2008, 2(2), 185–194.
  28. Lai, T. C., Morsi, Y. S., Das, S., & Owida, A. Numerical analysis of particle deposition in asymmetrical human upper airways under different inhalation cycles. Journal of Mechanics in Medicine and Biology 2013, 13(04), 1350068.
  29. Lee, D., Park, S. S., Ban-Weiss, G. A., Fanucchi, M. V., Plopper, C. G. & Wexler, A. S. Bifurcation model for characterization of pulmonary architecture. The Anatomical Record 2008, 291(4), 379–389.
  30. Lin, C., Lee, J.-H., & Hsieh, C.-M. The correlation between subcarinal angle and left atrial volume. Acta Cardiologica Sinica 2012, 28(4), 332–336.
  31. Liu, Y., So, R., & Zhang, C. Modeling the bifurcating flow in a human lung airway. Journal of biomechanics 2002, 35(4), 465–473.
  32. Martonen, T., Yang, Y., & Xue, Z. Effects of carinal ridge shapes on lung airstreams. Aerosol science and technology 1994, 21(2), 119–136.
  33. Mason, E. C., Wu, Z., McGhee, S., Markley, J., Koenigs, M., Onwuka, A., & Zhao, K. Computational fluid dynamic modeling reveals nonlinear airway stress during trachea development. The Journal of Pediatrics 2021, 238, 324–328.
  34. Minnich, D. J., & Mathisen, D. J. Anatomy of the trachea, carina, and bronchi. Thoracic surgery clinics 2007, 17.4, 571–585.
  35. Morita, K., Takeishi, N., Wada, S., & Hatakeyama, T. Computational fluid dynamics assessment of congenital tracheal stenosis. Pediatric Surgery International 2022, 38(12), 1769–1776.
  36. Murray, J., Brown, A., Anagnostou, E., and Senior, R. Widening of the tracheal bifurcation on chest radiographs: value as a sign of left atrial enlargement. AJR. American journal of roentgenology 1995, 164(5), 1089–1092.
  37. Nucci, Gianluca; SUKI, Béla; Lutchen, Kenneth. Modeling airflow-related shear stress during heterogeneous constriction and mechanical ventilation. Journal of Applied Physiology 2003, vol. 95, no 1, p. 348-356.
  38. Pedley, T., Schroter, R., and Sudlow, M. Energy losses and pressure drop in models of human airways. Respiration physiology, 1970; 9, 3, 371–386.
  39. Pedley, T., Schroter, R., and Sudlow, M. Flow and pressure drop in systems of repeatedly branching tubes. Flow and pressure drop in systems of repeatedly branching tubes. Journal of Fluid Mechanics 1971, 46(02), 365–383.
  40. Piemjaiswang, R., Shiratori, S., Chaiwatanarat, T., Piumsomboon, P., & Chalermsinsuwan, B. Computational fluid dynamics simulation of full breathing cycle for aerosol deposition in trachea: effect of breathing frequency. Journal of the Taiwan Institute of Chemical Engineers, 97, 66-79.
  41. Ren, S., Li, W., Wang, L., Shi, Y., Cai, M., Hao, L., & Luo, Z. Numerical analysis of airway mucus clearance effectiveness using assisted coughing techniques. Scientific reports 2020, 10(1), 1–10.
  42. Sahni, D., Batra, Y. K., & Rajeev, S. Anatomical dimensions of trachea, main bronchi, subcarinal and bronchial angles in fetuses measured ex vivo. Pediatric Anesthesia 2006 , 18(11), 1029–1034.
  43. Singhal, S., Henderson, R., Horsfield, K., Harding, K., & Cumming, G. Morphometry of the human pulmonary arterial tree. Circulation Research 1973, 33(2), 190–197.
  44. Schroter, R. and Sudlow, M. Flow patterns in models of the human bronchial airways. Respiration physiology 1969, 7(3), 341–355. [CrossRef]
  45. Tran, T. M., Huh, S., Kim, S., Cui, X., & Choi, S. Numerical investigation of the effect of tracheostomy on flow and particle transport characteristics in human airways. Physics of Fluids 2022, 34(12), 121901.
  46. Tsega, E. G., Katiyar, V. K., & Gupta, P. Numerical Simulation of Transport and Deposition of Dust Particles in Human Tracheobronchial Airways. International Journal of Biomedical Science and Engineering 2019, 7(1), 8.
  47. Wall, W. A., & Rabczuk, T. Fluid–structure interaction in lower airways of CT-based lung geometries. International Journal for Numerical Methods in Fluids 2008, 57(5), 653-675.
  48. Wang, Y., Liu, Y., Sun, X., Yu, S., & Gao, F. Numerical analysis of respiratory flow patterns within human upper airway. Acta Mechanica Sinica 2009, 25(6), 737-746.
  49. Wang, W., Dai, Z., Li, J., & Zhou, L. A hybrid Laplace transform finite analytic method for solving transport problems with large Peclet and Courant numbers. Computers & geosciences 2009, 49, 182–189.
  50. Wang, X.-Q., Mujumdar, A. S., & Yap, C. Effect of bifurcation angle in tree-shaped microchannel networks. Journal of Applied Physics 2007, 102(7), 073530.
  51. Weibel, E. R.; Cournand, A. F.; Richards, D. W. Morphometry of the human lung, 1st ed.; Springer Verlag OHG Berlin Goltigen Heidelberg: Berlin, Germany, 1963. [Google Scholar]
  52. West, J. B. Respiratory physiology: the essentials, 9th ed.; Lippincott Williams & Wilkins: United States, 2012. [Google Scholar]
  53. Liu, W., Sun, H., Lai, D., Xue, Y., Kabanshi, A., & Hu, S. Performance of fast fluid dynamics with a semi-Lagrangian scheme and an implicit upwind scheme in simulating indoor/outdoor airflow. Building and Environment 2022, 207, 108477.
  54. Van Ertbruggen, C., Hirsch, C., & Paiva, M. Anatomically based three-dimensional model of airways to simulate flow and particle transport using computational fluid dynamics. Journal of Applied Physiology 2005, 98(3), 970–980.
  55. Xia, G., Tawhai, M. H., Hoffman, E. A., & Lin, C.-L. Airway wall stiffening increases peak wall shear stress: a fluid-structure interaction study in rigid and compliant airways. Annals of biomedical engineering 2010, 38(5), 1836–1853.
  56. Xu, X., Wu, J., Weng, W., & Fu, M. Investigation of inhalation and exhalation flow pattern in a realistic human upper airway model by PIV experiments and CFD simulations. Biomechanics and modeling in mechanobiology 2020, 19(5), 1679–1695.
  57. Yang, X., Liu, Y., So, R., & Yang, J. (2006). The effect of inlet velocity profile on the bifurcation COPD airway flow. Computers in biology and medicine 2006, 36(2), 181–194.
  58. Yu, G., Zhang, Z., & Lessmann, R. Computer simulation of the flow field and particle deposition by diffusion in a 3-D human airway bifurcation. Aerosol Science and Technology 1996, 25(3), 338-352.
  59. Zhang, W., Xiang, Y., Lu, C., Ou, C., & Deng, Q. Numerical modeling of particle deposition in the conducting airways of asthmatic children. Medical Engineering & Physics 2020, 76, 40-46.
  60. Zhang, Z., Kleinstreuer, C., & Kim, C. S. Effects of curved inlet tubes on air flow and particle deposition in bifurcating lung models. Journal of Biomechanics 2001, 34(5), 659-669.
  61. Zhao, Y. and Lieber, B. B. Steady inspiratory flow in a model symmetric bifurcation. Transactions of the ASME-K-Journal of Biomechanical Engineering 1994, 116(4), 488–496.
  62. Zhao, Y., & Lieber, B. B. Steady expiratory flow in a model symmetric bifurcation. (1994).
  63. Zobaer, T., & Sutradhar, A. Modeling the effect of tumor compression on airflow dynamics in trachea using contact simulation and CFD analysis. Computers in Biology and Medicine 2021, 135, 104574.
Preprints 68251 g001
Figure 1. Geometric parameters of Weibel model [51].
Figure 1. Geometric parameters of Weibel model [51].
Preprints 68251 g001
Preprints 68251 g002
Figure 2. Bifurcation model geometry proposed by Lee et al [29].
Figure 2. Bifurcation model geometry proposed by Lee et al [29].
Preprints 68251 g002
Preprints 68251 g003
Figure 3. Computational grid for (a) BA’s study and (b) CRR’s study.
Figure 3. Computational grid for (a) BA’s study and (b) CRR’s study.
Preprints 68251 g003
Preprints 68251 g004
Figure 4. Grid detail near the branch wall.
Figure 4. Grid detail near the branch wall.
Preprints 68251 g004
Preprints 68251 g005
Figure 5. Comparison of the axial velocity profiles with Zhao & Lieber experiments for inhalation [61] and exhalation [62].
Figure 5. Comparison of the axial velocity profiles with Zhao & Lieber experiments for inhalation [61] and exhalation [62].
Preprints 68251 g005
Preprints 68251 g006
Figure 6. Indicative scheme of the zones and points where the results are reported.
Figure 6. Indicative scheme of the zones and points where the results are reported.
Preprints 68251 g006
Preprints 68251 g007
Figure 7. Velocity profiles on coronal plane.
Figure 7. Velocity profiles on coronal plane.
Preprints 68251 g007
Preprints 68251 g008
Figure 8. Velocity profiles on sagittal plane.
Figure 8. Velocity profiles on sagittal plane.
Preprints 68251 g008
Preprints 68251 g009
Figure 9. Velocity profiles at exhalation stage for TRA branch
Figure 9. Velocity profiles at exhalation stage for TRA branch
Preprints 68251 g009
Preprints 68251 g010
Figure 10. Comparison of coherent structures for different Reynolds numbers. Inhalation stage with B A = 55 . 9 ∘ .
Figure 10. Comparison of coherent structures for different Reynolds numbers. Inhalation stage with B A = 55 . 9 ∘ .
Preprints 68251 g010
Preprints 68251 g011
Figure 11. Comparison of coherent structures for different BA’s. Inhalation stage at Re = 2000 .
Figure 11. Comparison of coherent structures for different BA’s. Inhalation stage at Re = 2000 .
Preprints 68251 g011
Preprints 68251 g012
Figure 12. V-W profile and secondary flow patterns for inhalation stage.
Figure 12. V-W profile and secondary flow patterns for inhalation stage.
Preprints 68251 g012
Preprints 68251 g013
Figure 13. Comparison of coherent structures for different Reynolds numbers. Exhalation stage with B A = 55 . 9 ∘
Figure 13. Comparison of coherent structures for different Reynolds numbers. Exhalation stage with B A = 55 . 9 ∘
Preprints 68251 g013
Preprints 68251 g014
Figure 14. Comparison of coherent structures for different BA’s. Exhalation stage at Re = 2000
Figure 14. Comparison of coherent structures for different BA’s. Exhalation stage at Re = 2000
Preprints 68251 g014
Preprints 68251 g015
Figure 15. V-W profile and secondary flow patterns for exhalation stage.
Figure 15. V-W profile and secondary flow patterns for exhalation stage.
Preprints 68251 g015
Preprints 68251 g016
Figure 16. Pressure drops across the axial axes of the branches.
Figure 16. Pressure drops across the axial axes of the branches.
Preprints 68251 g016
Preprints 68251 g017
Figure 17. Pressure drops across the axial axes of the branches at exhalation stage.
Figure 17. Pressure drops across the axial axes of the branches at exhalation stage.
Preprints 68251 g017
Preprints 68251 g018
Figure 18. Wall Shear Stress on inner wall for inhalation.
Figure 18. Wall Shear Stress on inner wall for inhalation.
Preprints 68251 g018
Preprints 68251 g019
Figure 19. Wall Shear Stress on outer wall for exhalation.
Figure 19. Wall Shear Stress on outer wall for exhalation.
Preprints 68251 g019
Preprints 68251 g020
Figure 20. Velocity profiles on coronal plane at 20 % of the LMB length.
Figure 20. Velocity profiles on coronal plane at 20 % of the LMB length.
Preprints 68251 g020
Preprints 68251 g021
Figure 21. Velocity profiles on coronal plane at 50 % of the LMB length.
Figure 21. Velocity profiles on coronal plane at 50 % of the LMB length.
Preprints 68251 g021
Preprints 68251 g022
Figure 22. Velocity profiles on sagittal plane. Taken at 20 % of the LMB length.
Figure 22. Velocity profiles on sagittal plane. Taken at 20 % of the LMB length.
Preprints 68251 g022
Preprints 68251 g023
Figure 23. Velocity profiles on sagittal plane. Taken at 50 % of the LMB length.
Figure 23. Velocity profiles on sagittal plane. Taken at 50 % of the LMB length.
Preprints 68251 g023
Preprints 68251 g024
Figure 24. Comparison of coherent structures for different Reynolds numbers. Inhalation stage with r c = 0 .
Figure 24. Comparison of coherent structures for different Reynolds numbers. Inhalation stage with r c = 0 .
Preprints 68251 g024
Preprints 68251 g025
Figure 25. Comparison of coherent structures for different CRR’s. Inhalation stage at Re = 2000 .
Figure 25. Comparison of coherent structures for different CRR’s. Inhalation stage at Re = 2000 .
Preprints 68251 g025
Preprints 68251 g026
Figure 26. Vorticity profile on coronal plane at 20 % of the LMB length.
Figure 26. Vorticity profile on coronal plane at 20 % of the LMB length.
Preprints 68251 g026
Preprints 68251 g027
Figure 27. V-W profile at coronal plane and secondary flow patterns in LMB airway.
Figure 27. V-W profile at coronal plane and secondary flow patterns in LMB airway.
Preprints 68251 g027
Preprints 68251 g028
Figure 28. Pressure drop across the airways axial line.
Figure 28. Pressure drop across the airways axial line.
Preprints 68251 g028
Preprints 68251 g029
Figure 29. Wall Shear Stress on inner wall for inhalation.
Figure 29. Wall Shear Stress on inner wall for inhalation.
Preprints 68251 g029
Table 1. Weibel’s model parameters in the first three generations.
Table 1. Weibel’s model parameters in the first three generations.
AIRWAY D (mm) L (mm)
TRA 18 120
LMB / RMB 12 47.6
LUB / RUB / LLB / RIB 8.3 41.5
Table 2. Subcarinal angles reported by Haskin & Goodman [19].
Table 2. Subcarinal angles reported by Haskin & Goodman [19].
AGE [YEARS] ANGLE [°]
21 - 30 64.2
31 - 40 58.9
41 - 50 61.9
51 - 60 63.1
61 - + 55.9
All ages 60.8
Table 3. Model Parameters values for CRR study.
Table 3. Model Parameters values for CRR study.
Branch D (mm) Len. (mm) α (°) Rd (mm)
TRA 16 80 35 81.5
LMB / RMB 14 70
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated