Air Pollution Control Benefit and Cost Assessment of ICAPC2013 Paper ID: A-1-08

Air Pollution Control Benefit and Cost Assessment of ICAPC2013 Paper ID: A-1-08

Computational fluid dynamicsmodeling ofairflow and aerosol depositionin realistichuman upper airway

HU Gui-lin1, Lin Jiang1, FAN Jian-ren2+, Pan Deng2

1. School of Light Industry, Zhejiang University of Science and Technology, Hangzhou 310023, P. R. China;

2. State Key Laboratory of Clean Energy Utilization, ZhejiangUniversity, Hangzhou 310027, P. R. China

ABSTRACT

Based on the CT (Computerized Tomography) scanned images of the respiratory tract from a 19-years-old boy, a realistic three-dimensional geometric model from nasal cavity to the upper six-generation bronchial is rebuilt in this paper. In order to investigate airflow and particle deposition in this complexgeometric structure, RNG k-ε turbulence model was used to describe the airflow structureand the aerosol particleswere tracked and analyzed in the frame of Lagrange. The velocity field and pressure field of airflow, local deposition position of particles in the realistic human upper respiratory tractfor three airflow rates such as 15L/min, 30L/min and60L/min were computed and discussed. The obtained complex airflow structure in model is important for the study of respiration principles. Results indicate that the airflow structure and pressure distribution is strongly affected by airflow rates and respiratory states as inhalation and exhalation. The same characteristic both for inhalation and exhalation is that the maximum velocity is located at pharynx part. Furthermore, the particle deposition position was visualized, which can provide reference for the accurate drug delivery. The regional deposition number of different diameters under different airflow rates was counted. It indicates that deposition in the nasal cavity occupy the major part of total deposition in the whole model. In the bronchia, the particle preferred deposition in the carinal ridge and downstream at the inner tube of the daughter tubes for the effect of impact between particles and wall. And the deposition in the bronchia proportionally increases as airflow rate.

KEYWORDS

airflow; aerosol deposition; realistic human upper airway; computational fluid dynamics modeling

INTRODUCTION

Particles ranging from toxic particulate matter to drug aerosols may be inhaled into people’s respiratory tract, and some of them can be either harmful or therapeutic to humans depending upon the aerosol material, deposition site and local concentration. In turn, these aspects and parameters are greatly determined by the airflow field, particleproperties, breathing pattern and geometric airway characteristics. Thus, good knowledge of airflow structure, deposition of particles in the upper respiratory tractcan provide credible theoretical foundation not only for preventing illness of respiratory tract effectively. On the other hand, it is desirable todeliver pharmaceutical aerosols to the targeted part accurately, and thuseffectively reduce drug dosage and pain of patients (Schlesinger and Lippmann, 1978; Balásházy et al., 2003).

Hitherto, especially in the last decade, many research literatures about computation fluid dynamics modeling of airflow and particles movement in the respiratory tract are available all over the world. Most models employ low-Reynolds k-ω model and standard k-ε model to describe gas flow, track particles movement in the frame of Lagrange. Using a mouth-to-trachea replica plus three generations of theWeibelAbronchia, Zhang et al. (2005) as well as Zhang et al. (2006) investigated numerically steady inhaled particle transport and deposition. Employing aWeibel-type bifurcation and considering chronic obstructive pulmonary disease (COPD),Yang et al. (2006) and Liu et al. (2003) compared the resulting airflow structures under different inlet conditions. Tang et al. (2004)presents a brief summary of flow features in the human respiratory system and simulates an airflow field based on a 3D real-anatomical geometry of the human nasal cavity. The dispersed phase in the flow field is solved by Lagrangian particle-tracking approach. Hofmann et al. (2001) computed the inspiratory deposition efficiencies of ultrafine particles with 6500 nm in airway Generations 3 and 4 for different inlet flow rates.In light of the importance of targeted drug-aerosol delivery, Li et al. (2007) study the relation between particle-release positions at the trachea inlet and particle depositions at specific lung sites, which is found to be greatly influenced by the complex airway geometry and the flow-rate magnitude. Just recently, the LES method is employed to the study of airflow in the human respiratory tract (Jin et al., 2007). For the geometric model, it can be found that most simulation is based on the idealized, i.e., planar and symmetric, lung airway model published by Weibel (1963). Recently, modern imaging techniques allowed for even more detailed mapping of the human respiratory system (Cebral and Summer, 2004; Sera et al., 2003). CFD simulations based on the CT scanned images are conducted (Xi and Longest, 2007; Shi et al., 2007).Croce et al. (2006) measured pressure–flow relationships in human plastinated specimen of both nasal cavities and maxillary sinuses, and compared them to those obtained by numerical airflow simulations in a numerical three-dimensional reconstruction issued from CT scans of the plastinated specimen.however, the bronchial structure is very limited for the difficulty of image resolution. Especially, the structureof nasal cavity is very complex; the rebuilt of it is relatively difficult.Nowak et al.(2003) conducted computational fluid dynamics (CFD) simulations of airflow and particle deposition in geometries representing thehuman tracheobronchial tree. Two geometries were used in this work; one of them is based on a CT scan of a cadaver lung cast. Flow conditions used included both steady-state inhalation and exhalation conditions as well as time-dependent breathing cycles. Particle trajectories were calculated in each of these models by solving the equations of motion of the particle for the deterministic portion of particle displacement. The trapping of particles on the wall surfaces was monitored, and the locations of trapping in each generation were recorded.

The objectives of this paper is to rebuilt a three-dimensional geometric model of realistic human upper respiratory tract includingnasal cavity, pharynx, larynx, trachea andthe uppersix-generation bronchial based on the CT scanned images.On the other hand, RNG k-ε turbulence model is employed to describe the airflow and track the particle in the frame of Langrange.And then the airflow structure and inhalable particles deposition position in this realistic human upper respiratory tract is obtained.

2 DESCRIPTION OF MODEL

2.1 GEOMETRIC MODEL

2.1.13D CT IMAGES SCANNING

The helical 3-D CT scanning was performed on a SOMATOM Sensation 16 station with the following parameters: 120 kV, 140 mAs, rotation time of 0.5 seconds, thickness of 0.75 mm and multislice scanning with 16 slices per rotation. Transverse CT imageswere obtained from a healthy 19-year-old adult male (the volunteer was awake and reclining). These images were taken through the full extent of respiratory tract spanning a vertical distance of approximately 50 cm while the volunteer was instructed to hold his breath under full inhalation condition in order to avoid false image between inhalation and exhalation.Thus, 690 images have been scanned for this research.

2.1.2. MODEL RECONSTRUCTION

The subsequent images data were then transferred to a computer workstation for analysis by image and graphics recognition technologies. With the CT data provided, we created a 3-dimensional geometric model of the upper airway. This 3-dimensional structure provides measurements of the upper airway at any specified spatial position. In the CT image, density values are represented by gray scale values. In our method, we identify a region of interest by defining a range of gray values. For the upper airway, the gray value is the lowest because that the air is filled in it. Using image recognition techniques as MIP and VRT, the upper airway including nasal cavity, pharynx, larynx, trachea and bronchia is converted into an STL model, which indicates that a triangle mesh will be generated around the selected volume. The obtained 3-dimension upper airway is shown in figure 1. The geometric structure of the upper airway is very complex, especially for the nasal cavity, which is different from most present model used by computational fluid dynamics. This will cause complicated flow fields into the upper airway, and complicated particle transport and deposition principles along the airway. Thus it is important base for further research. The more real flow characteristics can be obtained using the more real geometrical model.

Figure.1 The three-dimensional geometrical model of realistic human upper respiratory tract rebuilt from the CT scanned images including nasal cavity, pharynx, larynx, trachea and upper six generation bronchia

The trachea and bronchia are separately illustrated in figure 4 for more particular knowledge of the structure. From which we can find that the bronchia is asymmetric; especially the center lines of the bronchia form a spatial solid angle.

2.2 GOVERNINGEQUATION

The airflow structure will be very complex because sharp shrink and expand surface exist in the pharynx-larynx, bifurcate structure exist in the bronchia and the nasal cavity with more complex structure.Thus, the RNG turbulent model is adopted in this paper, which provides an option to account for the effects of swirl or rotation by modifying the turbulent viscosity appropriately. Furthermore the RNG model is more responsive to the effects of rapid strain and streamlines curvature than the standard k- ε model. The continuity and movement equations:

[1]

[2]

The turbulence kinetic energy, k, and its rate of dissipation, ε, are obtained from the following transport equations:

[3]

[4]

[5]

Where i、j＝1、2、3, ρ is fluid density, t is time, is fluid velocity vector, x denotes the spatial coordinate, stands for average pressure, is fluid viscosity, is turbulent viscosity, is the gravity in the i direction, kis turbulence kinetic energy, εis turbulence dissipation rate, represents the generation of turbulence kinetic energy due to the mean velocity gradients,The quantities and are the inverse effective Prandtl numbers for kandε,βis thermal expansion coefficient, Eij implies main time-average strain rate.andare model constants. The model constants and have values derived analytically by the RNG theory, used by default in this paper, are

,

The particle phase in the air can be considered as dilute phase, and neglecting the interaction between particles. The Lagrangian dispersed phase model is used for the prediction of the trajectory of a particle. This is done by integrating the force balance on the particle. In a Lagrangian reference frame, movement vector of a single particle can be written as:

[6]

Where, mp is particle mass, dp is particle diameter, virepresents particle velocity vector, CDstands for particle resistance coefficient, g acceleration of gravity. The three terms on the right hand are Stokes resistance, gravity force and Saffman’s lift force considered in this paper.

Particle velocity can be obtained by integrated equation (4), furthermore particle position can be tracked by following equation:

[7]

The particle deposit efficiency (DE) is introduced to quantificationally analyze particle deposit in different region of the respiratory tract:

[8]

Where Nddenotes the number of particle deposited on wall of someregion in the respiratory tract, Ntrepresents the total number of particle inhaled into the model.

2.3 INITIALCONDITIONANDBOUNDARYCONDITION

The pressure outlet boundary is imposed on outlet of bronchial, i.e., the pressure is zero and gradient of other variables is zero. When the particles impact with model wall,it is marked as deposition and the track computation is over, which accords with the real instance well.Different from many papers, the pressure inlet condition is selected in this paper. One box extending from nasal inlet is added to consider the inlet effect, and the pressure at nasal inlet is not equal one atmosphere but the boundary of the box.

2.4 Numerical procedures

In this paper, block partition grid method is employed in order to save computer source and computation time. Thus, most regions of this complex computational domain are partitioned by tetrahedron grids as figure 2, and the total number of grid is not very large. In order to exclude effect of the grid number on computational results, the grid number about 200,000, 250,000, and 300,000 are tested.It is found that 250,000 grids is enough accurate for this problem resolution, and adopted for this paper.

Figure2.Mesh of the computational domain

The numerical solutions of the continuity and momentum equations as well as the k and ε equations were carried out with Fluent 6.2.In the present simulation, the coupled airflow equationsand pressure fields are solved by SIMPLEC (Semi-Implicit Method for Pressure-Linked Equations Consistent) algorithm with under-relaxation. The particle transport equation was also solved with the commercial software in combination with particle tracking. Particles were released simultaneously at the box boundary after the constant flow field has been established. The force of particle caused by airflow is considered, which is called one-way method.

3. RESULTSANDDISCUSSIONS

3.1 AIRFLOWSTRUCTURES

Three airflow rateswith Q＝15L/min (low-level breathing at resting condition), Q=30L/min(middle-level breathing atlight exercise condition), Q=60L/min (high-level breathing atmoderateexercise condition) are adopted as computational case respectively. The airflow fieldat several typical across sections is discussed in the following section.

3.1.1AIRFLOWFIELDDURINGINHALATIONPHASE

Figure3 illustrates the velocity contours in the pharynx-larynx part and bronchia bifurcation for breath intensity Q=60L/min during inhalation phase. Be similar to cases Q=15 and 30 L/min, the maximum velocity of the whole modelwith Vmax=3.82 m/s still lies in the pharynx. Similarly, there has one shrink structure in the larynx; the maximum velocity at the cross-section in the larynxexhibits close to the posterior wall. After entering into bronchia, the maximum velocity at the cross-section in the bronchiagradually moves towards the center of the tube, and the gradient of pressure and velocity is small. Totally, the velocity in the whole airway is two times higher than that of breath intensity Q=30L/min. The velocity of bifurcation is plotted in figure 14b. It can be seen from figure that the gas stream erodes the carinal ridge more intensively as the increase of gas velocity. Thus, the deposition of the gas stream carrying particles will be larger than that of breath intensity Q=30L/min. Similarly, the carnial ridge of every generation bronchia and its daughter tube will be the hot spot for particles deposition.

(a)Velocity contours in the pharynx-larynx / (b)Velocity contours in the bronchia bifurcation
Figure. 3. Velocity contours in the pharynx-larynx and bronchia bifurcation for breath intensity Q=60L/min during inhalation phase

3.1.2AIRFLOWFIELDDURINGEXHALATIONPHASE

Figure4depicts the velocity contours in the pharynx-larynx and bronchia bifurcation for breath intensity Q=15L/min during exhalation. The velocity contours in the pharynx-larynx is shown figure17a, from which it can be seen that there have relatively higher velocity in the contraction segment of pharynx and larynx, and the maximum velocity with 3.84m/sin the whole model is located at the pharynx part, which is the same as inhalation phase. Different from inhalation phase, there has no vortex close to the posterior wall of larynx for the geometry structure of the joint between the trachea and larynx.The maximum velocity at the cross-section in the contraction segment doesn’t lie at the posterior wall. The gradient of velocity close to the center position is small because the human bronchial is smooth, and then the gas stream become stable gradually during exhalation. The velocity in the trachea bifurcation is illustrated in figure 17b.It can be seen from figure that there have counteract and superpose process for the partial velocity with same direction at the joint of same-generation two bronchial. After fully uniting, the gas stream attainsmaximum velocity at the center of trachea.

(a)Velocity contours in the pharynx-larynx / (b)Velocity contours in the bronchia bifurcation
Figure 4. Velocity contours in the pharynx-larynx and bronchia bifurcation for breath intensity Q=15L/min during exhalation phase

3.2 PARTICLESDEPOSITION

In order to investigate the particle deposition position and hot spots, Lagrange method is implied to track the particle movement. The particle deposition at different parts and under different breath intensities are visualized in figure5. The figure a plots the deposition in the bronchia, from which we can see that the particle preferred deposition in the carinal ridge and downstream at the inner tube of the daughter tubes for the effect of impact between particle and wall. The deposition in the nasal cavity, pharynx-larynx, and trachea is shown in figure. It can be seen from figure that the deposition in the front of the pharynx-larynx is small for the air flow affected by the geometry structure.

(a) Nasal cavity (b) pharynx-larynx (c) trachea

(d) Bronchial (e) bifurcation

Figure. 5. Distribution of particles deposition in different part of the respiratory tract with d=1μm and Q=15 L/min during inhalation phase, the darkish particle indicates deposition at the rear wall

4. CONCLUSIONS

A three-dimension geometric model of a real human upper respiratory tract including nasal cavity, pharynx, larynx, trachea and upper six generations bronchial based on the CT scanned images is rebuilt in this paper. RNG k-ε turbulence model is used to compute airflow, particle movement and deposition efficiency in the obtained complex structure geometrical model. The three-dimensional flow structureunder three different breathing intensity15L/min, 30 L/min and 60 L/min are computed and analyzed. Then particle of 1μm diameter in every part of respiratory tract are tracked and visualized under Lagrange frame. The result indicates that the total deposition efficiency in the whole model slackly increases as the increase of breathing intensity. The deposition efficiency in the nasal cavity decreases as the breath intensity increases; however, the deposition efficiency in the pharynx-larynx, trachea and bronchia slackly increases as the breath intensity increases. Among all the parts, the deposition efficiency in the nasal cavity is the largest, and the deposition efficiency in the trachea is the smallest.

REFERENCES