Influence of Morphological Parameters on the Flow Development Within Human Airways

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}^{\circ }$ and ${65}^{\circ }$for all ages is described by Haskin & Goodman [19], meanwhile in the study carried out by Bipinchandra et al [5] a range between ${50}^{\circ }$ to ${130}^{\circ }$ was found. Furthermore, Christou et al [9] reports ranges from $65.{04}^{\circ }$ to $122.{01}^{\circ }$ for men and from $69.{46}^{\circ }$ to $113.{94}^{\circ }$ for women. On the other hand, the study by Sahni et al [42] showed that this parameter has values ranging from ${42}^{\circ }$ to ${75}^{\circ }$ 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^{\circ }$ and $63.1^{\circ }$ 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}^{\circ }$ 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 $\alpha$ describe the curvature of the bifurcation. The CRR was non-dimensionalised using the following expression,

$$r_c=\frac{CRR}{r_t}$$

(1)

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,

$$\nabla ·\overrightarrow{\mathit{v}}=0$$

(2)

$$\frac{\partial {\mathit{u}}_i}{\partial t}+\nabla ·\left({\mathit{u}}_i\overrightarrow{\mathit{v}}\right)=-\frac{1}{\rho }\nabla p+\nabla ·(\nu \nabla {\mathit{u}}_i)$$

(3)

where ${\mathit{u}}_{\mathit{i}}$ and $\overrightarrow{v}$ are the i-th velocity component and the velocity vector field, respectively Equally, in Equation (3), $\mathit{p}$ is the fluid pressure and $\nu =\mu /\rho$ its kinematic viscosity with $\rho$ standing for the density of the fluid. For the purposes of this work, the Reynolds number has been defined as

$$\mathit{Re}=\frac{U_{ave}D}{\nu }$$

(4)

Where $U_{ave}=0.5U_{max}$. Air properties were defined at ambient temperature of 15^o, with $\rho =1.23[kg/m^3]$ and $\nu =1.48e-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\left(r\right)=U_{max}\left(1-{\left(\frac{r}{R}\right)}^2\right)$$

(5)

Where $U\left(r\right)$ is the inlet velocity, $U_{max}$ is the maximun velocity, r a radial position and R the inlet radius. Human breathing rates range from $0.2L/s$ to $2.5L/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 $\mathit{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 $\mathit{Re}=500$ and $\mathit{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

$$\frac{\partial \overrightarrow{v}}{\partial n}=0$$

(6)

In the surface of the wall a Dirichlet no-slip boundary condition was applied to the velocity field at all airway walls,

$$\overrightarrow{v}=0$$

(7)

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\times {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\times {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\times {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^{+}=\frac{u_{\tau }y}{\nu }$$

(8)

where $\nu$ is the fluid’s kinematic viscosity, $u_{\tau }$ is the friction velocity defined as

$$u_{\tau }=\sqrt{\frac{{\tau }_w}{\rho }}$$

(9)

the variable y is the wall-normal distance from the cell centre to the bounding wall, and ${\tau }_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 ${\tau }_w$ was calculated at the maximun velocity of our test set, i.e. $\mathit{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_{max}$ for $\mathit{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 $\mathit{Re}=500$ and $\mathit{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 $\mathit{Re}$ numbers. For the $\mathit{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 $\mathit{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 $\mathit{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 $\mathit{Re}=500$ and $\mathit{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 $\mathit{Re}=500$, and Figure 8d and Figure 8f for $\mathit{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 $\mathit{Re}=500$ and are more noticeable for lower BA values. The velocity profiles in the coronal and sagittal planes for $\mathit{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

$$\mathit{Q}=\frac{1}{2}{\left(\right||\mathbf{\Omega }||}^2-{\left|\right|\mathit{S}\left|\right|}^2)$$

(10)

where $\mathbf{\Omega }$ is the rotation tensor and S is the rate of strain tensor. The dimensionless factor $Q^{*}=Q\left(D^2\right)/{\left(U_{bulk}\right)}^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_{bulk}$. The two vortex cores generated after each bifurcation point for all branches are shown in Figure 10a for $\mathit{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 $\mathit{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 $BA=55.9^{\circ }$ are similar to those for $BA={70}^{\circ }$ at $\mathit{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^{\circ }$. 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 $\mathit{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 $\mathit{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 $\mathit{Re}=2000$. The comparison of these structures for $BA=55.9^{\circ }$ and $BA={70}^{\circ }$ 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 $\mathit{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 $\mathit{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 $\mathit{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 $\Delta 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 $BA=55.9^{\circ }$ and $BA={70}^{\circ }$ is up to 20% for the first bifurcation point.

The pressure drop profiles in the exhalation state for $\mathit{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 $BA=55.9^{\circ }$. The decrease of the WSS between the angle of $BA=55.9^{\circ }$ and $BA={70}^{\circ }$ 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 $BA=63.1^{\circ }$ 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 $\mathit{Re}=500$ and $\mathit{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 $\mathit{Re}=500$ and $\mathit{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 $\mathit{Re}=500$ and $\mathit{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 $\mathit{Re}=500$ and $\mathit{Re}=2000$. These curves were built using a normalized length d/D and a normalized vorticity $\omega /{\omega }_{max}$ 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 $\mathit{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 $\mathit{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 $\mathit{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 $\mathit{Re}=500$ and $\mathit{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 $\mathit{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 $\mathit{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 $\mathit{Re}=500$ to $\mathit{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.