Pressure Loss Coefficients for Asymmetric Bifurcations of Pulmonary Airways with Predetermined Flow Distributions

Computational Fluid Dynamics simulations of inspiratory airflow in asymmetric bifurcations have been performed in order to determine the influence of the asymmetry and Reynolds number on pressure losses over a physiologically relevant range for pulmonary airways; thus the results of this work can contribute to the understanding of respiratory ventilation in health and in disease. A key a priori insight to the design of the study is that the flow distribution in respiratory bifurcations can be largely independent of the local losses; and therefore, is predetermined by the boundary conditions in these calculations. The results, presented in the form of pressure loss coefficients, indicate that asymmetry and downstream conditions are significant for severe restrictions and laminar flow; but are relatively insignificant for turbulent flow conditions and for flow through the healthy branch. Pressure Loss Coefficients for Asymmetric Bifurcations of Pulmonary Airways with Predetermined Flow Distributions


Introduction
The distribution of inhaled air throughout human lungs is influenced in part by the pressure losses in lung airways. These losses are reflections of patient health; alterations in measured pressure drops can indicate the presence of respiratory disease [1]. The seminal work of Olson and colleagues [2] clearly demonstrated the contributions of various anatomical features to these energy losses, isolating the effects of branching angle and cross-sectional area on calculated pressure losses in pulmonary airways.
Pressure losses also have implications for the effective treatment of respiratory diseases, as they affect the extent to which inhaled medication is distributed throughout the lung. It is therefore of interest to examine the pressure losses within diseased airways. Heterogeneous airway obstructions are associated with a variety of respiratory diseases, including asthma, COPD, and cystic fibrosis, and can limit ventilation of obstructed regions of the lung [3][4][5][6].
Due to the complexity of multi-generational branching in the lungs, it is difficult to develop a single model that accurately and reliably predicts the flow of inhaled gas in all circumstances. Following the single-patient model of Weibel [7], researchers developing geometric models have had to choose whether to model one patient's airways accurately, or to simplify the morphology in order to generalize their results. Lung airway shape and structure vary significantly both within and among subjects [8,9], the average branching angle varies from 15° to 70° [10], and the length-to-diameter ratio depends strongly on both generation and parent diameter [11]. Most previous airflow models have simplified the morphology of the lungs by assuming that each bifurcation is symmetric [12]; others have included more complexity, but have not modeled airway obstructions [11,13].
It is the intention of the current work to use Computational Fluid Dynamics (CFD) to demonstrate and quantify the influence of asymmetry and bifurcation geometry on the pressure losses in branching airways, while recognizing that the ultimate flow distribution is largely determined by the particular placement of the bifurcation in the complex network of airways within the lung. For example, the flow distribution is often much more influenced by downstream lung volume or heterogeneous lung compliance than the asymmetry of a particular branch. Thus, a parametric analysis is performed where the morphologies and flow conditions are varied independently. Airway diameters are varied in such a way as to simulate the effect of natural asymmetry or the reduced lumen that can occur due to some respiratory diseases. The results in the form of pressure loss coefficients are provided and general observations are discussed that could inform the development of future airway fluid mechanics models.

Numerical method
The airflow through the lungs is governed by the Navier-Stokes equations, and the flow in any airway can be characterized by the nondimensional Reynolds number, Re: Where ρ and μ represent the air's density and viscosity, respectively; v is the mean velocity, and D is the diameter of the airway of interest. The air is assumed to have constant properties: density ρ=1.225 kg/m 3 and viscosity μ=1.7894*10-5 kg/s m. In the current work, the reference Reynolds number is calculated using conditions at the inlet of the parent airway. changes direction from inspiration to exhalation, the current work considers only inspiratory flow, and examines steady flow at a range of Re values to examine the full range of physiological flow conditions. The three-dimensional bifurcation geometries of interest were meshed using GAMBIT (ANSYS, PA, USA), and the flow was simulated using the finite volume solver FLUENT (ANSYS, PA, USA). A mesh refinement study confirmed that the pressure drops and loss coefficients were sufficiently well resolved. A "boundary layer" mesh, six cells deep, was used at the walls to better quantify shear stresses, and to achieve grid independence of the solutions. The number of nodes and elements in the mesh varied depending on the diameter and length of the daughter vessels; the average number of elements for all models was 835,000. Within FLUENT, second-order-accurate discretization schemes were used for all terms. Pressure-velocity coupling was achieved using the SIMPLE algorithm, and the transitional k-kl-ω model, a three-equation eddy-viscosity model for laminar and turbulent kinetic energies (k and kl, respectively) as well as inverse turbulent time scale (ω), was incorporated. The model is based on two transport equations, one for intermittency and one for the transition onset criteria in terms of momentum thickness Reynolds number. The transport equations are intended for the implementation of correlation-based models into general-purpose CFD methods [14]. The theoretical framework and correlation constants for the turbulence model are beyond the scope of this paper but are provided in the FLUENT Theory Guide [15]. The Reynolds-Averaged Navier-Stokes (RANS) model devolves to the laminar Navier-Stokes equations when the turbulent fluctuations are negligible such that it was applicable for all of the cases considered, even those that remained laminar throughout.
The boundary conditions included no-slip and no-penetration at the walls; a prescribed uniform inlet velocity for the parent airway; and outlet conditions on the daughter vessels to simulate the influence of downstream effects in the pulmonary network. The outlet boundaries were given classic outflow conditions on downstream derivatives but included flow rate weighting to prescribe the distribution of inlet flow to the daughter vessels, such that each daughter received a specified portion of the parent vessel's volumetric airflow.

Calculation of losses
In addition to the mass and momentum conservation equations that govern flow, it is possible to consider the energy balance for steady, incompressible flow from any location A to location B along the airway: Where p and v are the average cross-section pressure and velocity at each station, and z is the position with respect to gravity g. The kinetic energy term is modified by a coefficient that varies to account for different flow profiles: α=1 for a blunt velocity profile and α=2 for parabolic. The term on the right hand side of equation (2), known as the head loss, represents the sum of all the resistive energy losses per unit mass as the flow moves from A to B. The energy balance in equation (2) assumes that there are no other sources of energy (heat transfer, or moving walls) present.
The head loss can also be expressed as the combination of major losses (for the straight parent and relevant daughter segments) and minor losses K (for the bifurcation): Here f is a friction factor whose value depends on the flow regime and Reynolds number Re. For simplicity we will refer to the loss coefficient K without the subscript because it is the only one considered herein. For laminar, fully developed flow, f is analytically determined: For turbulent flow, f is calculated from the Blasius correlation: It should be noted that the energy balance of equation (2) is applied from location A along the parent vessel to a station B on either daughter vessel, so that there are two values of head loss associated with a single bifurcation. These values will be equal if the bifurcation is completely symmetric, but in the case of asymmetric bifurcation, there will be two values of h A-B and, accordingly, two values of K. We will use the notation "left" and "right" to distinguish between the loss coefficients for the two daughter vessels.
The desired coefficient K for either daughter vessel can be isolated by rearranging the previous two expressions.

Model geometry
To explore flow patterns over the full range of flow rates and inhaled gases for many generations of airways of the human tracheobronchial system would require an immense effort; herein a single base model was employed where the inlet Re was varied from Re=100 (which occurs in lower regions of the lung) to Re=2000 (which occurs nearer the trachea) [16]. A schematic is shown in Figure 1. For the base model the diameter of the parent vessel was held at 0.01 m and the non-restricted right daughter at 0.008 m for all simulations, while the diameter of the left daughter branch was varied to explore the effects of downstream constriction. It should be emphasized that while FLUENT requires the input of particular dimensional data, the reformulated inputs and outputs as discussed below are applicable throughout the lung. The remaining geometric parameters were held constant. For example, the bifurcation angle is fixed at 70°, with each daughter vessel having a branching angle of 35° from the parent vessel's mid plane and the inlet and outlet lengths where 20D to allow for flow development .
As indicated above, the effects of disease on flow in the lungs are modeled by flow boundary conditions and by changes to the geometric model. Three key non-dimensional parameters are varied independently: the Reynolds number Re, defined in equation (1); and the ratios of daughter to parent volumetric flow rates (equation 7) and diameters (equation 8). The volumetric flow rate ratio Q ratio is expressed: Where Q daughter represents the left daughter branch's proportion of the inlet flow; the flow rate in the right daughter vessel is prescribed to be 1-Q daughter . The value of Q ratio is permitted to vary from 0.1 to 1 to simulate a range of disease states. The ratio of daughter to parent vessel diameters is: Where D daughter represents the left daughter vessel, as the right daughter vessel diameter is held fixed (at 0.8D parent ) throughout the study. D ratio is varied from 0.4 to 1.0. The entry and exit terms of Equation 6 were taken at one diameter upstream (L/D=-1) and 19 diameters downstream (L/D=19) of the bifurcation to allow for fully developed flow. In all, 100 CFD simulations were performed (Table  1) with two K values being calculated per simulation, one for the left (diseased) branch and one for the right (healthy) branch.

Results and Discussion
The calculated flow patterns are consistent with previous studies of flow in bifurcations [13,[17][18][19]. Flow separation occurs in some cases, with the imposed flow distribution (Q ratio ) being the strongest predictor of separation. This is seen in Figures 2 and 3: in Figure 2, a pathological Q ratio of 0.8 results in significant separation in the right daughter vessel; in Figure 3, a Q ratio of 0.3 yields low pressure but no backflow or separation in the same region.
The pressure loss coefficients in the form of K values for the left (diseased) and right (healthy) branches are given in Table 1 in the Supplementary Material. Some general comments can be made. The left (diseased) branch K values increase significantly when flow is shunted from the right (healthy) branch. For example, for D ratio =0.4 (a 50% reduction in diameter from the base case) and Re=100, the K value increases over 400% (from 5.2 to 28.2) as the Q ratio increases from 0.5 (symmetric) to 1.0 (totally shunted). However, this is the case only for laminar flow. Consider the equivalent case (D ratio =0.4, Q ratio increases from 0.5 to 1.0) with Re=2000; here the K value actually decreases from 2.4 to 1.0 (but note that the overall pressure drop still increases). Another example where the presence of turbulence appears to dampen the increase in K value can be observed by comparing simulations 18-20 where the parent Re=2000. K values increase only about 40% (from 7 to 10.1). But note that there is a relative maximum because for Q ratio =0.8 K=16.1. This can be explained by the fact that for the daughter Re=1993 and 2492 (i.e., the Re straddles the turbulent limit) for Q ratio =0.8 and 1.0, respectively; i.e., the presence of turbulence results in a relatively lower and more uniform value independent of the other parameters.
Over the full parameter range examined the K values for the right (healthy) daughter are relatively uniform (average of 1.3, standard deviation of 0.4) and thus close to the symmetric cases (Simulation Nos. 53, 58, 63, 68, 73). A similar observation as discussed above can be made for the left (diseased) daughter for R=2000 (average of 1.2, standard deviation of 0.5). In fact, even for the laminar cases (Re=100, 200, 500 and 1000) when the flow is shunted away from the diseased to the healthy branch (Q ratio =0.1, 0.3 and 1.0) the K values average 1.6, though the standard deviation is somewhat greater, 4.7. These results suggest that for a wide range of cases the effects of downstream conditions and asymmetry are relatively negligible and the symmetric values can most likely be employed giving reasonable results. However, for more severe restrictions in laminar flow downstream conditions should be taken into account.
We attempted without success to find analytical correlations for these results to make them easier to apply directly into models of airway pressure loss. However, given the observations stated above, perhaps limited correlations for only laminar flow and severe disease may be found that have adequate accuracy.
The CFD simulations performed in the current work reflect physiological conditions, and do not prescribe fully developed flow. The resulting pressure drops are compared with those calculated using fully developed flow models for three representative cases in Figure  4. Although fully developed assumptions are implicit in the pressure drops and loss coefficients generally used in large network models, the differences shown in Figure 4 indicate that for shorter daughter sections that typically occur in the lung the inertial loss due to the bifurcation  will be somewhat overestimated. For example, for the cases illustrated in Panels A, B and C of Figure 4, if the daughter branch ended at 5 diameters instead of 19 as indicated the errors would be 23%, 9%, and 30%, respectively.

Conclusion
Pressure loss coefficients have been provided for asymmetric lung bifurcations over a physiological relevant range of incoming Reynolds number, the diameter ratio of parent to daughter, and the flow ratio. The flow ratio was prescribed to account for the fact that individual bifurcations are part of complex lung airway networks.