The influence of boundary conditions to the flow through model of upper part of human respiratory system

Respiratory system represents relatively large system of gradually branching channels which can be hardly solved by numerical simulations. Nowadays, research in this area is focused to solve problems in selected parts of respiratory tract rather than whole system. This simplification comes with problem of accurate assessment of boundary conditions on model geometry. Geometry used on Department of Thermomechanics and Environmental Engineering on Brno University of Technology consists of mouth cavity, larynx, trachea and bronchial tree up to seventh generation of branching. This article is focused on comparison of two different settings of boundary conditions steady inspiration during light activity regime. First set of boundary conditions represents commonly used setting with zero pressure resistance on outlet from the model and second method deals with more realistic assumption, where incomplete 3D geometry is coupled with the rest of bronchial tree described by 1D equations and also correlated by the amount of air, which flows in specific lung lobe. The article observed differences in individual mass flow through the model branches under different conditions and its influence on the flow structures.


Introduction
For the human being, breathing is vital process which secures intake of ambient air into alveolar part of the lungs where blood stream is oxygenated and carbon dioxide is discharged outside the human body.During this phenomena, the air is fed through oral or nasal cavity to pharynx, where continues to trachea and to the lungs.The lungs are paired asymmetrical organ which represents the system of gradually branching channels who ends in alveolar sacs.These sacs are further connected with bloodstream [1].Inspiration then leads to enlargement of lung volume through the movement of diaphragm and intercostal muscles and due to different pressure between air inside the lungs and air outside the human body.Negative pressure in the lungs causes airflow inside the body.The airways in the bronchial tree are grouped in five separate areas called lung lobes, which may be determined using techniques of visualisation of the respiratory tract by the computed tomography or magnetic resonance [2].Estimation of the size of these lobes is important, because the amount of inhaled air is divided by these values and values of the volume change of the lobes during breathing.If we want to properly set the boundary conditions of the numerical simulation of transport and deposition in human upper airways, it is appropriate to take account of the above facts.Since the respiratory tract is extensive and it would be computationally demanding to solve whole of it from mouth to the alveolar sacs, it is usual to solve only portion of the system that is interesting for the specific research [3,4].However, it is difficult to prescribe a realistic boundary condition at the parts which was interrupted due to reduction of the size of the respiratory system.Typical solution of this situation includes the prescription of the pressure boundary conditions on the outlets from the interrupted parts of the bronchi region of the model.Determination of the correct value of pressure resistance which would represent a real pressure loss of the following rest of bronchial tree is difficult.In the past, these losses were not considered and a pressure boundary conditions were prescribed with zero pressure resistance [5].Way to provide more realistic boundary conditions is to couple the 3D solution with much easier 1D model of the downstream of bronchial tree, based on idealized geometry of the lung ( [6,7] for example).This method was developed by Fry et al [8] and nowadays evolved by Grandmont et al [9].Their approach is based on Weibel symmetric model and pressure losses solved by Hagen-Poiseuille law equation.The Weibel model does not consider the segmentation of the lungs to the lung lobes which was taken into account during our research.

Model
Geometry of the human upper airways (see figure 1) containing the oral cavity assembled with a suction nozzle of cylindrical shape, larynx, trachea and bronchial tree up to the seventh generation of branching of conductive paths of bronchial tree was used for the model.This geometry consists of several different parts obtained as a result of cooperation with worldwide research centers.Segment of oral cavity, larynx and part of trachea was obtained from Lovelace Respiratory Research Institute in Albuquerque, USA [10] as a wax cast and then digitalized using 3D scanner ATOS into stereolithographic format.Part of the trachea was taken by computed tomography at the University Hospital St. Anna in Brno and model of the bronchial tree was obtained in MeVis computer data format from the Institute of Anatomy and Cell Biology at the Justus-Liebig University in Giessen, Germany [11].The geometry of the bronchial tree ends in the seventh generation of branching.Segmentation of the models bronchial tree is described on the

Numerical setup
The article deals with two cases of the calculation.First option deals with usual prescription of boundary conditions, where the inlet to the model was prescribed by velocity boundary condition which corresponds to the respiratory condition (light activity which means 30 l/min through the inlet of diameter of 20 mm which corresponds to the speed of 1.59 m/sec), and at the outlets of the bronchial tree the pressure boundary condition with zero pressure resistance were prescribed.This approach was already used in past researches [12] and is usual in cases where calculations of transport and deposition of aerosols in respiratory tract were solved.Second option is based on the knowledge that the lungs are divided into five independent units (lung lobes) and airflow through the branches is caused by pressure differences between inner and outer environment of the lobes.Amount of the air that reaches the lungs is given by the breathing regimes, which depends on the man activity and this air is divided into lungs according to the size of the lung lobes.Sizes of the lobes were obtained from the review of articles focused on measuring of the lung lobes volume using computed tomography or magnetic resonance.Thanks to this information, coefficient equal to the amount of the air which flows to the specific lobe in dependence of the total amount of inhaled air, was assigned to these lobes (see table 1.).Then the values from specific subjects were averaged for each lobe and used for calculation.Thanks to the knowledge of the location of the boundary conditions (see Figure 1.) in the bronchial tree model, we know in which lung lobe the specific boundary condition lead and we will be able to divide inhaled amount of the air by these coefficients.Next step is based on the distribution of inhaled air to the specific boundary condition in a particular lobe.Certain number of boundary conditions can be found in the specific lobe, where velocity inlet condition is prescribed.Because the branches of interrupted parts of bronchial tree, where the inlet boundary conditions can be found, are located between the fifth and the seventh generation of branching of bronchial tree, the correlation based on pressure resistance of the following part of bronchial tree had to be done.The required pressure resistance was found by coupling of 1D equation of pressure loss calculation for laminar flow in the cylindrical pipe to the CFD calculation of flow in presented model, where Poiseuille law (1) depending on velocity at the outlet from the model was calculated for the rest of bronchial tree represented by Wiebel model geometry.
Differences in pressure losses was used as a correlation factor, which was later prescribed as a final boundary condition on outlets from the model geometry.

Calculation
The calculation was performed using the CD-Adapco StarCCM+ version 8.02.008.On the geometry, the polyhedral mesh with density of 3.5 million cell containing prismatic layer of eight rows for better description of boundary layer, was created.Numerical simulation was made using uRANS with the K-omega SST (Menter) model of turbulence with low Reynolds number treatment activated.This approach was selected according to work of Kleinstreuer et al. [17].Selected points of interest were monitored during calculation with respect to monitor the oscillation of velocity values.The calculation was considered as solved, when convergence of residual falls below 10 -4 .

Results
Mass flow of air through selected parts of the model was monitored during calculation.By these monitoring, mass flow values on constricted branches of the models bronchial tree was determined.These values were then counted on the basis of mass conservation law and air distribution into each lung lobe was determined.Particular values of mass flow into specific lung lobe can be found in table 2.  2 shows uneven distribution of air into each lung lobes for both variants of the calculation.Graphic comparison of the values from table 2 and data obtained from literature review is shown in figure 2. In the case of variant A, the lowest distribution of air can be seen in the upper left lobe while the largest amount of mass flow can be found in lower right lung lobe.If we assume that the inhalation of the air is driven from the models inlet nozzle in mouth cavity, the airflow is influenced by the local resistance of individual parts of the bronchial tree.Bronchial tree is a set of gradually branching channels which contains a large number of bifurcations, where incoming air is divided onto two separate air streams that goes to different parts of the bronchial tree.Each bifurcation represents embedded local pressure resistance, and greater amount of air flows to the path where pressure resistance is lowest.For example, the first bifurcation provides distribution of air flow into left or right part of the lungs.In case where inhalation is initiated from mouth nozzle, the branch which leads into the right lung puts less pressure resistance to the air flowing through it thanks to its higher diameter ( right branch has 13 mm, left branch has 10 mm ) and is also diverted from the axis of the trachea by a lower angle ( right = 48°, left = 68°).Higher amount of air flows into right lung than to left lung thanks to this parameters.Then the air is further distributed into each lung lobes by the pressure resistance of the branches that forms the lung lobe.At least the air flowing into the right middle lobe, which is the smallest of all lung lobes.While in the case of Option A, flows into the left lung 30% of the inhaled air, while the second variant is 44.9%, which leads to a uniform distribution of air within the entire lung.The ratios between the mass flow of air which flows into lung lobes for the left or right lung are similar while use of both variants.

Figure 3 Velocity profile in trachea
Figure 3 shows a comparison of velocity profiles in line probe passing through the axis of trachea, 2 cm above carinal ridge of the first bifurcation, where the division of the inhaled air into two streams leading to the left and to the right lung can be found.Comparison shows that the method has not influence on the shape of the velocity profile, only a change in the ratio between the speeds of the left and right side of the trachea is visible.The air velocity at the right side of the trachea (left side of the graph), where the pressure resistance on entry to the right lung was higher during the calculation using the method A, than when calculating method B.

Conclusion
Realistic distribution of flow is suitable for future calculations of aerosol transport and deposition as the amount of air that flows into the individual branches of the bronchial tree which affects the distribution of velocity profile and the number of particles who enters into specific parts of the lung.Calculation above shows that the different approach to the prescription of boundary conditions based on coupling with 1D equation and using of correlations as information about size of each lung lobes is suitable to ensure a realistic distribution of air to the lungs.

Figure 1 .
For the purposes of the research, the model is divided on five segments which represents airway paths into different lobar lobe.These lobes are five, two on the left (sinister) side and three on the right (dexter) side of the model with respect of view of the subject.Name of the parts are follows: Right upper lobe (RUL), right middle lobe (RML), right lower lobe (RLL), left upper lobe (LUL) and left lower lobe (LLL).

Figure 2
Figure 2 Comparison between calculated variants

Table 1 .
Coefficients of lung lobe sizes

Table 1 .
Mass flow into specific lung lobe