Numerical investigation of airflow, heat transfer and particle deposition for oral breathing in a realistic human upper airway model

Author’s Accepted Manuscript Numerical investigation of airflow, heat transfer and particle deposition for oral breathing in a realistic human upper airway model X.Y. Xu, S.J. Ni, M. Fu, X. Zheng, N. Luo, W.G. Weng www.elsevier.com/locate/jtherbio

PII: DOI: Reference:

S0306-4565(16)30438-7 http://dx.doi.org/10.1016/j.jtherbio.2017.05.003 TB1935

To appear in: Journal of Thermal Biology Received date: 31 December 2016 Revised date: 30 April 2017 Accepted date: 17 May 2017 Cite this article as: X.Y. Xu, S.J. Ni, M. Fu, X. Zheng, N. Luo and W.G. Weng, Numerical investigation of airflow, heat transfer and particle deposition for oral breathing in a realistic human upper airway model, Journal of Thermal Biology, http://dx.doi.org/10.1016/j.jtherbio.2017.05.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Numerical investigation of airflow, heat transfer and particle deposition for oral breathing in a realistic human upper airway model X.Y. Xua, S.J. Nia, M. Fub, X. Zhenga, N. Luoa, W.G. Wenga, * a

Institute of Public Safety Research, Department of Engineering Physics, Tsinghua University, Beijing, 100084, China

b

Hefei Institute for Public Safety Research, Tsinghua University, Hefei, Anhui Province, 320601, China

*The corresponding author: Wenguo Weng Institute of Public Safety Research, Department of Engineering Physics, Tsinghua University, Beijing, 100084, China. E-mail: [email protected] Tel: +86-10-62796917 Fax: +86-10-62792863

1

Abstract Inhalation injury from exposure to fire smoke is one of the causes of burn-related death. In this study, a realistic three-dimensional human upper airway model was built from magnetic resonance imaging (MRI) scanned images, including the nasal, oral, pharynx, larynx, trachea and part of the first generation of the tracheobronchial tree, as well as a tissue region from the pharynx to the upper bronchi. The Transition Shear Stress Transport (SST-transition) turbulence model, Pennes bioheat transfer equation, convective boundary conditions and a Lagrangian frame were applied and verified with experimental measurements to simulate the airflow fields, temperature distributions and particle deposition in the human airway model. The effects of flow rate, inhalation temperature and particle diameter were studied. It showed that the oral cavity is more likely to be affected by the inlet air conditions. The mucosa in the oral, pharynx and larynx are more likely to cause the thermal injury. The inspiration flow rate significantly influences the airflow fields, temperature distributions and particle deposition fraction interior of the human upper airway model, especially in the pharynx-larynx region. The rising flow rate, inhalation air temperature and particle diameter all contribute to boosting the total deposition fraction in the model. The heated particles with a higher temperature are more likely to be deposited in the oral cavity and the influence of the inlet temperature has a minor influence in the case of a bigger particle diameter. Keywords: thermal burn; human upper airway; inhalation injury; CFD; airflow; particle deposition

1. Introduction Inhalation injury is one of common effects of a fire smoke accident. It is a leading cause of burning morbidity and mortality, by directly damaging the human upper airway through heat, the lower airway through chemical injury and/or the inhalation of toxic gas at the site of the fire (You et al., 2014, Watt et al., 2016). Furthermore, the smoke-related inhalation injury is often connected to noxious chemicals, and the toxic injury usually occurs in the trachea and lung tissues (Desai et al., 1989). 2

Computational fluid dynamics (CFD) simulations for drug delivery and toxic particles deposition on the airflow and particle deposition in different parts of human respiratory tract have been studied (Kleinstreuer and Zhang, 2010, Park and Wexler, 2008, Wang et al., 2009). Fluid flow and deposition were simulated in a replica of the human mouth-throat employing the laminar, standard k-ε turbulence model, detached eddy simulation (DES) and large eddy simulation (LES) methods, and concluded that for the simulation of medication aerosols inhaled at a steady flow, LES and DES provided more accurate results than the RANS turbulence model (Stapleton et al., 2000, Jayaraju et al., 2008). Low-Reynolds number (LRN) SST k-w turbulence model simulation was conducted and found to be consistent with the experiments conducted by Phase Doppler Anemometry (PDA) in almost all the cross sections (Elcner et al., 2016). LES and three widely used RANS turbulence models were compared to evaluate the flow fields in a locally constricted conduit and an idealized mouth-throat model, with the study of nanoparticle deposition for the accurate lung deposition concentrations of therapeutic nanoparticles, assuming that the SST transition model gave a better prediction of turbulence kinetic energy profiles in some cases while LES required more computational time than RANS turbulence models in spite of providing instantaneous velocity fluctuations (Zhang and Kleinstreuer, 2003b, Zhang and Kleinstreuer, 2011a). Micron-particle transport with realistic inlet conditions was studied in a (WeibelTypeA) 16-generation model (Zhang et al., 2009) and nasal airway (Inthavong et al., 2012) to determine the inhalation toxicology effects of exposure to atmospheric particles. Micron and nano-size particle deposition were simulated in the human lung with a semi-empirical model of single breath transport and deposition to estimate particle dosimetry in humans to evaluate long-term PM exposure (Park and Wexler, 2008). Nanoparticles deposition was shown to be influenced by the breathing route (from nasal to oral breathing), the outlet flow rate ratio and the mass transfer coefficients in a nasal-oral-tracheobronchial airway (Zhang and Kleinstreuer, 2011b). Fibrous particles, wood dust and particles with various morphologies were also analyzed in a realistic human nasal cavity to understand the deposition pattern (Tian et al., 2007, Inthavong et al., 2008, Inthavong et al., 2009a). Besides, studies of the heat transfer in the human respiratory tract are limited. A simplified human upper airway model focused on the deposition fractions of nano-scale particles, temperature distributions and velocity fields inhaling cold air and normal air with the thermal boundary condition of a constant wall temperature (Zhang and Kleinstreuer, 2003a). Then 3

a fluid and heat transfer on the human nasal cavity was conducted to consider the indoor humidity under cold and normal air conditions (Inthavong et al., 2009b). An averaged convective heat transfer coefficient distribution of the human airway under steady and transient breathing was calculated in indoor environments via inhalation in two different models and summarized as functions of breathing airflow rates (Phuong et al., 2016). A simulation of velocity and temperature distributions of inhalation thermal injury was generated in a human upper airway and provided the correlations for the Nusselt number for each individual part of the airway, which regarded the wall temperature as the constant experimental figure (Goodarzi-Ardakani et al., 2016). A more realistic model considering the individual’s physiological status, such as the metabolic rate, thermal conductive, water evaporation and blood perfusion rate, was constructed to evaluate the injury caused by air during the early stage of a fire (Lv et al., 2006, Xu et al., 2015). In vivo and in vitro experimental measurements were also carried out to validate the simulation models and results. Studies of the human upper airway flow in the vitro measurements were conducted, using the Particle Image Velocimetry (PIV) measurement method to accurately measure air flow fields studies (Kobayashi et al., 2010). A comparison of the PIV data and the CFD results of a nasal-oral upper human upper airway model had achieved relatively good agreement in the trachea region (Phuong and Ito, 2015). Endoscopic PIV was used to measure the sagittal plane velocity fields of an idealized model of the human oropharynx under steady inspiration (Heenan et al., 2003). Anatomical measurements (Lin et al., 2001), a metered dose inhaler (Hsu and Swift, 1999), fluorescence spectrometry (Cheng et al., 1999), gamma scintigraphy (Grgic et al., 2004) and the advanced positron emission tomography (PET) (Lizal et al., 2015) methods have also been applied to measure diverse particles in the human mouth and throat deposition. Since there are some restrictions on using humans for the in vivo heat transfer experiments, canine models are mostly used as substitutes. The airway surface temperatures measured in most experiments are considered to be correlated with thermal injury severity (Zhao et al., 2013a). An experiment on canine models showed that the air temperature decreases markedly from 80-320 ℃ at the oral cavity to 39.3-137.9 oCat the lower margin of the vocal fold, causing a more severe injury to the larynx and the lower airway. (Rong et al., 2011). Blood circulation plays an important role in heat absorption with the increase of the regional blood flow 4

rate, which leaves the airway tissue less injured than the skin (Wan et al., 2016, Zhao et al., 2013a). However, thermal injury to the respiratory tract is usually limited to dry air and is rarely due to combined burn and smoke inhalation injury. Particles of the smoke, such as the micron carbon particles, which are mainly affected by the airflow field, carry lots of heat during their transportation with deposition and movement in the human respiratory tract. Therefore, particle deposition can significantly affect the accuracy of the thermal injury evaluation. In addition, the realistic respiratory tract geometric model and the idealized ones can be quite different, and the shape of the model can greatly affect the flow fields. It’s quite important to use a realistic human upper airway model, precisely study the deposition position, deposition fraction of the fire smoke and the temperature distributions for thermal injury evaluation. In this study, a three-dimensional geometric human upper airway model, from nasal, oral to some part of the first generation of the tracheobronchial tree, was built from the magnetic resonance imaging (MRI) scanned images of a heathy female. Steady airflow fields conducted by the SST-transition turbulence model, heat transfer of human physiological status and particle deposition analysis in the frame of Lagrange were performed and validated. The effects of respiratory flow rate, inhalation temperature and particle size on airflow fields, temperature distributions and particle deposition fractions were discussed.

2. Methods 2.1. Airway geometry The original respiratory tract data were developed from the MRI scans of a healthy, 23-year-old, nonsmoking Asian female volunteer, with a body mass index (BMI) of approximately 22. The MRI scans produced 451 cross sections of the airway. Digital Imaging and Communications in Medicine (DICOM) were used to transfer and store the medical images. The subsequent MRI images were transferred to a 3-dimensional structure of the airway by using Amira® (VSG), to generate and modify 3D surface models from 2D contour data. After the process of segmenting, modifying, smoothing and merging, the geometric structure of the airway was generated as a STL 5

file, as shown in Fig. 1(a). It was composed of the nasal and oral cavity, pharynx, larynx and G0 (trachea)-part of G1. The airway model has a height and volume of 24.83 cm and 82.97 cm3, respectively. An elliptic cylinder which fitted the shape of the larynx and trachea region was generated to be the simplified tissue part around the respiratory tract, as shown in Fig. 1(b). The thickness was around 3cm. As the structure of the nasal and oral cavity part is complex, the simplified tissue part was not added.

Fig. 1. (a) Schematic diagram of the human upper airway model; (b) Schematic diagram of the human upper airway model and simplified tissue 2.2. Modelling of Airflow and heat transfer The airflow structure in the human upper airway is quite complicated as the range of the human breathing rate is wide from 10L per minute (LPM) calm rest to 90 LPM under strenuous exercise. It may show up the laminar, laminar-to-turbulent transition and turbulent airflow. With the gradual decrease in the cross-sectional area of the mouth-throat, nasopharynx and pharynx-larynx, the velocity of the airflow increases rapidly and causes a transition from the laminar to the turbulence, and in some areas may even result in separation of the boundary layer. To accurately reflect this laminar-to-turbulence transfer in a numerical simulation, the SST transitional model was used for the respiratory flow (Zhang and Kleinstreuer, 2011a, Inthavong et al., 2012). For reasons of brevity, the turbulence model equations for the SST-transitional model (Menter et al., 2006) are 6

not reproduced here. For the transportation of heat in the respiratory tract, the model below the larynx region can be simplified as 2 parts, the air layer and the tissue around the wall, shown in Fig. 2. The thickness of the tissue is around 3 cm, and the theoretical model can be derived with the conservation of energy.

Fig. 2. Schematic diagram for the simplified upper airway For the air temperature, the energy equation can be written as

  k g Tg (  g uiTg )  ( ) xi xi cg xi

(1)

where  g , cg , k g , Tg are the density, specific heat, thermal conductivity and the calculated temperature of the air, respectively. As the temperature difference affects the buoyancy force, the Boussinesq model is employed for the buoyancy force.

f j   g [1  T (Tg  Tref )]

(2)

where T , Tref are the thermal expansion coefficient and buoyancy reference temperature, respectively. g =9.81m / s

2

is the gravitational acceleration. 7

For the tissue area, based on the Pennes bioheat transfer equation, the resulting equation is as follows

T  (kt t )  Wb cb (Ta  Tt )  Qm  0 xi xi

(3)

where kt is the tissue thermal conductivity; cb denotes the specific heat of blood; Wb the blood perfusion; Qm the metabolic heat generation rate; Tt the tissue temperature and

Ta =37℃ the arterial blood temperature is assumed as a constant. The effect of metabolic heat generation and blood circulation are considered as source terms and implemented by a user-defined function (UDF). During the thermal process, at the tissue–air interface ((5), (11)), the generalized boundary condition (BC) for the heat transfer is convection. The thermal boundary condition for the air-wall interface ((1), (6), (7), (12)) was simplified, assuming that the airway tissue can provide or extract all the heat due to its rich blood supply, i.e., the wall exchanging heat with a constant surrounding temperature by convection. To solve the problem, the heat flux was assumed to approach zero at the deep tissue ((3), (9)) and at both ends of the tissue ((2), (4), (8), (10)),which is realistic for a biological body. The thermoregulation mechanisms of the biological bodies have been ignored because of the slight temperature increase induced. Three ambient conditions were considered in this paper. The inlet air temperature Tin  45C,80C,120C , which indicates the early stage temperature of a fire. At the outlet the backflow temperature was set at 37 oC. The boundary conditions can be written as follows (5), (11): kt

T =h f (Tt  Tg ) y

(1), (6), (7), (12): k g

(3), (9):

Tg y

(4a)

=h f ' (Tg  T0 )

(4b)

Tt =0 y

(2), (4), (8), (10):

(4c)

Tt =0 x

(4d) 8

where h f is the heat convection coefficient between the tissue and the flowing air stream,

T0 =37℃ is the assumed surrounding temperature and h f ' is the heat convection coefficient between the air wall and the surrounding. 2.3. Particle equations The Lagrangian particle tracking method was used for the dispersion of the particle trajectory, because the particle phase is generally a dilute phase of the two-phase flow in the respiratory tract. The particle is assumed to be spherical with a uniform radius. The particle dispersions were calculated to either reach a boundary surface and be trapped or escape from the outlet. Ignoring the collision between the particles, the movement of the particles can be determined by Newton's second law of motion:

mp

du p dt

 FD  m p

g ( p  g )

p

 Fs

(5) where m p , u p and  p are the mass, velocity and density of the particle, respectively. u g and

 g is the velocity and density of the air, respectively. FD represents the drag force and Fs is the other possible forces, such as the Saffman force, the virtual mass force, pressure gradient force, Basset force and thermophoretic force. The drag force can be given as

FD 

3 g CD Re p 4  p d p2

(6) where d p is the particle volume equivalent diameter and  g is the molecular viscosity of the air. The particle Reynolds number, Re p and the drag coefficient CD are defined as

Re p =

 p d p u p  ug g

(7)

9

CD =a1 

a a2  32 Re p Re p

(8) where Equation (8) is the drag model for the spherical model, depended on Re p (Morsi and Alexander, 1972). As for the fire smoke, micron particles where

 p   g (Inthavong et al., 2012), the virtual

mass force, pressure gradient force, and Basset force can be ignored. But the Saffman force and thermophoretic force should be considered in the boundary area. Fs can be given as

1 2

Fs =1.6   g  g  d p ug  u p

dug dx pi

1 2



6 d p  g 2Cs ( K '  Ct Kn)

1 Tg (9)  g (1  3Cm Kn)(1  2 K  2Ct Kn) m pTg x pi '

where Ct  2.18 , Cm  1.14 are the momentum exchange coefficient and Cs  1.17 is the thermal slip coefficient. x pi is the particle position coordinate. Kn 

number, as  is the mean air free path. K  '

kg kp

2 is the Knudsen dp

, where k p is the thermal conductivity of the

particle. The Discrete Random Walk (DRW) model was also used, as the motion of the particles were affected by the random discrete turbulent eddies in the turbulence dispersion. The instantaneous velocity of the fluid consists of two parts, the average speed ui where i  x, y, z ,

ug  ux2  u y2  uz2 and the fluctuation velocity is ui ' . For the isotropic area, ui ' can be calculated by (Gosman and Loannides, 1983)

ui ' =

2 k 3

(10) where  is a random number subject to a normal distribution, k is the turbulent kinetic 10

energy. Due to the assumption of turbulence isotropy, the fluctuating velocities calculated by DRW can be higher than the actual values, which may then overpredict the particle deposition in some cases (Matida, et al., 2004). A near-wall correction was used to simulate the near-wall particle '

trajectories, with the fluctuating velocity normal to the wall un being expressed as

un ' =f v

2 k 3

(11) with

f v  1  e0.02 y



(12) where f v is a damping function component normal to the wall considering the anisotropy of

turbulence near the wall, which is used for y   60 (Zhang et al., 2005, Inthavong et al., 2012). The regional deposition of the micron particle in the human upper airways can be quantified in the terms of the deposition fraction ( DFp ) or deposition efficiency in a certain region ( DE p ), which can be defined as

DFp 

Nd 100% Nt

(13a)

DE p 

Nd 100% Nr

(13b) where N d is the number of trapped particles in a specific region , N t is the number of particles entering the human upper airway and N r is the number of particles entering this region. 2.4. Numerical methods 11

The numerical solution of the governing equations was carried out with Ansys-Fluent v17.2. The SIMPLE algorithm was used to solve the flow equations, and the second-order upwind scheme was employed to improve the accuracy. The mouth was set as the uniform steady velocity inlet and the nostrils were the wall boundaries with a wall temperature the same as the inspiratory temperature. The outlets were set as zero pressure outlets. The wall boundaries were all set as no-slip. For oral inhalation, the Reynolds number was 3592, based on the inlet velocity which was 5.03m/s for the steady inhalation (30 LPM). All the inlet Reynolds numbers showed that the airflow was not the turbulent condition. The effect of metabolic heat generation and blood circulation were considered as source terms and added by a user-defined function (UDF) into the tissue part. The variables used are summarized in Table 1 (Lv et al., 2006, Phuong et al., 2016).

Table 1 The used variables used in simulation. Variable

Value

Cb

4000 J/kg·℃

 tissue

1000 kg/m3

Qm

420 W/m3

K

0.5 W/m·℃

Ta

37 ℃

Wb

0.5kg/m3s

hf

30 W/m2·℃

hf '

50 W/m2·℃

Cg

1005 J/kg·℃

Kg

0.025 W/m·℃ 12

T

0.03676

The computational mesh was generated with Ansys ICEM-CFD. A tetrahedral unstructured mesh filled the airway passage and the tissue region. Prism layers were applied to the bounding respiratory walls, which contain the viscous sub-layers and can deal with any geometric features there, as shown in Fig. 3. The first grid point above the whole surface wall was given a value of 

the wall units ( y ), less than or equal to 1.0, and this criterion is maintained for all the simulations. The final mesh contains about 8.6 million and 1 million cells for the human upper airway model and the tissue part, respectively.

Fig. 3 Representative mesh outline of numerical airway model and tissue part The visualization and generation of the mesh, and the simulations were all conducted on a workstation with 128 Gb RAM, 2.4GHz CPUs and 16 processor cores. The solution of the flow field was assumed to have converged when the monitored variables, such as the dimensionless mass residual ( (total mass residual)/(mass flow rate)<10-4) were negligible.

3. Model validations To accurately simulate the fluid-particle dynamics, the validations from the experimental results of the airflow, particle deposition and inhaling air heat distribution are necessary. 3.1. Airflow 13

The simulated results of the airflow are validated with the measured non-dimensional scalar velocity distributions in the case of the 10L/min flow, given in Fig. 4. y  y /

A where A is the

cross-sectional area and ubulk represents the mean velocity in this cross section. There are some discrepancies in the geometry of the airways between the current simulation and the experiment data it is compared with (Johnstone et al., 2004). The flow velocity in the human upper airway changes considerably with the structure of the airway. Therefore, the results are comparable only in the trachea zone in the current validation, where the diameter is almost uniform and the airflow approaches the laminar. It can be seen that the simulated results of the SST transition model agree well with experimental data in terms of values and trend, better than that of the RNG k-ε model. The small discrepancies between the simulation and the experimental measurement can be caused by the differences of the geometry and the compared location, as well as by the stationary hot-wire probes used in the experimental measurement which doesn’t work well with the reverse flow (Ball et al., 2008).

Fig. 4. Comparison of simulated velocity profiles at trachea with the experimental data (Johnstone et al., 2004) at Q=10L/min. 3.2. Particle deposition 14

Fig. 5 provides comparisons of the micron particle deposition efficiency between the simulation results in the present upper airway model with oral inhalation and the in vitro experimental data of several different oral airway models (Grgic et al., 2004). In the mouth and throat deposition, the inlet conditions greatly affect the deposition efficiency. With the different geometries and inlet diameter conditions of the models, the Stokes number is employed to take the relevant length and velocity into account, calculated as

 d p 2U St  18 D

(14)

where D is the corresponding length scale for the special geometry or inlet conditions, U is the velocity scale, and  is the air dynamic viscosity. Compared with the experiments’ models, the bronchia part G1 of the present model is not included and the deposition efficiency of the simulation results doesn’t contain this part. It can be seen that the numerical prediction of the micron particle deposition is in agreement with the experimental results in different kinds of human oral airway models. The present discrepancies might be concerned with the different shape and the different inlet velocity direction of the present model.

15

Fig. 5. Comparison of the simulated particle deposition efficiency with the in vitro experimental data (Grgic et al., 2004). S1a, S2, S3, S4, S5 are marked as the experimental conditions.

3.3. Temperature distributions of inhaled hot air

Fig. 6. Comparison of predicted air temperature in trachea with the experimental measurement (Rong et al., 2011, Zhao et al., 2013b).

16

Fig. 6 compares the simulation results of the average temperature distributions in the trachea with the measurement of the inhalational thermal injury in the canine model (Rong et al., 2011, Zhao et al., 2013b). A periodic sine function of oral inhalation and exhalation for 4min and breathing the indoor air (20℃) for 10s every 1min, which is similar to the experiment, was conducted in the simulation. The breathing cycle was 4s and the inlet velocity was loaded as UDF. It can be seen that the simulated results are in reasonable agreement with the measured data. The compared discrepancies could be attributed to: (i) the experiment considered the blood temperature increase with the inhaled hot air temperature while it is as a constant in the simulation; (ii) evaporation exists in the mucous layer in the experimental canines and the equations are included in the computing air-tissue interface; and (iii) the device used to measure the air temperature in the trachea only gave a single point result, which might cause some experimental deviations. In the case of inhaling the hot air under 160℃, which may not reach the inhalation air hyperthermia but can represent the early stage of a fire, the blood temperature will not rise more than 1.5℃ (White and Cabanac, 1995, Sandsund et al., 2007, Zhao et al., 2013b) and it is reasonable for it to be assumed as a constant. In addition, circulational heat dissipation can be assumed to play an important role in the heat-absorbing process of the upper airway when inhaled air is less than 160 o

C and the mucous evaporation can be ignored.

In summary, the reasonably good agreements between the computational predictions and experimental results show that the present model can accurately analyze the laminar to turbulent fluid flow, particle deposition and heat transfer in the human upper airway for inhalation injury assessment.

4. Effects of respiratory flow rate, inhalation temperature and particles size Simulation results of airflow, particle deposition and temperature distributions are discussed for different inspiratory flow rates, inlet air temperatures and particle sizes. The particles diameters for fire smoke are bigger than 1μm, and less than 100 μm. In the current study, three kinds of particles with diameter (3, 5 and 6.5 μm), three inlet air temperatures (45 oC, 80 oC and 120 oC) and three airflow rates (10, 30 and 60L/min) are studied for the particle deposition in three regions—oral cavity, pharynx-larynx and trachea. 17

4.1. Airflow fields Fig. 7 shows the coded streamlines as well as dimensionless velocity contours at selected cross sections of the human upper airway model for oral breathing only under 30L/min flow rate and 80℃ inlet temperature, i.e., medium-level breathing in a heated air situation. There is a speed reduction after the air inhalation, and the velocity beneath the oral cavity is higher than that at the top. From the oropharynx to the pharynx/larynx region (from slice AA’ to DD’), the elevated airflow from the anterior part moves towards the low-speed region at the posterior of the human upper airway model, which enhances the airstream expansion and mixing in the pharynx/larynx region. With the narrowing of the cross-section area of the vocal fold, a jet appears in the laryngeal region, reaching a 2.9 times higher maximum velocity than that of the entrance. As the influence of the glottis nozzle structure, the high-speed zone shifts to the anterior side from the posterior, where a reverse flow occurs in the posterior. After the airstream flows into the trachea, it turns into a uniform and laminar flow with the Reynolds number at about 3500. The flow patterns provide a potential deposition place for the inhaled particles. As the particles are transported by the airflow, the spots of the peak velocities and changing flow directions suggest a high potential of particle deposition, such as the entrance of the oral cavity and the posterior of the oropharynx, where the flow direction changes nearly 50o.

Fig. 7. 3-D view of streamline contours as well as normalized velocity contours at selected cross sections in the upper airway model under the flow rate of 30L/min oral inhalation for the inlet air 18

temperature of 80 oC. The flow features for different inspiratory flow rates resemble those for Qin  30 L / min . However, the inspiratory flow rates can influence the detailed flow fields in the human upper airway model. Fig. 8 displays the normalized scalar velocity profiles at the central sagittal plane of six selected cross-sectional locations (as those from slice AA’ to FF’ in Fig. 7) under three oral breathing conditions, 10 L/min (light breathing), 30 L/min (medial breathing) and 60 L/min (heavy breathing). Although the breathing flow rate does not change much of the flow patterns for the oral inhalation, it affects the flow features near the wall in those selected cross sections. The higher the inspiratory flow rate, the more complexity velocity profiles appear in the oropharynx part, where the airstream goes to the posterior and a larger vortex appears. The turbulent fluctuations close to the wall become stronger after the constriction of the glottis, and then slightly weaker in the trachea. It can also be observed that the flow pattern of 30, 60L/min flows are more uniform than that in the case of the 10L/min flow.

Fig. 8. Comparisons of normalized velocity profile in the selected cross-sections of different breathing flow rates for the inlet temperature of 80 oC with only oral breathing. Fig. 9 shows the comparison of the normalized velocity profile in the selected cross-section of different inlet temperatures (45 oC, 80 oC and 120 oC) under the breathing flow rate of 10 L/min 19

oral breathing. It can be seen that the buoyant flow only slightly enhances the air mixing near the wall, and the inlet temperatures have a stronger impact on the oropharynx part in contrast with those on the trachea, where the blood circulation in the tissue removes much of the heat and reduces the temperature differences. It illustrates that the different inlet temperatures tend to have only minor influences on the airflow fields compared to the inspiratory flow rate. However, the thermal effects tend to be minor for the high flow rate cases, such as 30L/min and 60L/min flows, which are not shown in figures.

Fig. 9. Comparison of normalized velocity profile in AA’ cross-section of different inlet temperatures (45 oC, 80 oC and 120 oC) under the breathing flow rate of 10L/min oral breathing. 4.2. Temperature distributions Fig. 10 presents the temperature distributions in terms of non-dimensional temperature profiles, i.e.,  =(T  Ta )/(Tin  Ta ) in the central sagittal plane and different cross-sectional locations (as those from slice AA’ to FF’ in Fig. 7) in the human upper airway model under the 60L/min inspiration flow and the 120 oC -inlet temperature. It is shown that the heat spreads slightly to the nasopharynx with the fluxion of the airstream and the air temperature in the human upper airway 20

gradually decreases as the distance away from the entrance increases. The temperature in the tissue shows that the air-tissue interface temperature is above the other parts. The surface temperatures of the larynx and trachea region become quite uniform and much lower with the blood circulation effects in the tissue. Clearly, the temperatures of the posterior, where observed with the relatively high-velocity, are higher than those of the interior in the larynx (from slice BB’ to DD’), indicating the delaying of the blood heat absorption and the difficult heat dissipating with a high velocity. Meanwhile, in the trachea regions, the posterior air-tissue interface temperature is slightly lower than that of the interior and the heat absorption becomes more remarkable in the lower trachea, as shown in slice GG’. It is suggested that the posterior surface of the larynx and the interior interface areas of the upper trachea might more easily incur thermal injury.

Fig. 10. Non-dimensional temperature contours at selected cross-sectional slices under the flow rate of 60L/min oral breathing conditions for the inlet temperature of 120 oC. Fig. 11 displays the non-dimensional temperature profiles of six selected cross-sections (same as those in Fig. 8) under different oral breathing conditions, i.e. 10, 30 and 60L/min. Comparing the temperature profiles of Fig. 11(b-f) with Fig. 11(a), the convection with tissue, involving blood circulation effects, promotes more uniform temperature distributions under diverse flow rates. The higher the respiratory flow rate, the faster the air warms up, for a higher velocity implies a short residence time for the heat exchange between the air, tissue and blood. It can be expected that the 21

inspiration flow rate might cause the heat transfer in the human upper airway, leading to the thermal injury.

Fig. 11. Comparisons of non-dimensional temperature profile in the selected cross-sections of different breathing flow rate under oral breathing conditions for the inlet temperature of 45 oC.

Fig. 12. Comparison of non-dimensional temperature profile in the selected cross-section of different inlet temperatures (45 oC, 80 oC and 120 oC) under the breathing flow rate of 10L/min oral breathing conditions. 22

However, in the near wall areas surrounded with tissue, the dimensionless parameter θ is observed to slightly change with variations in inlet air temperatures (45 oC, 80 oC and 120 oC) under the same oral inhalation rate, seen in Fig. 12(b-f). It appears that the θ increases with a lower inlet air temperature, and the discrepancy of θ near the wall has the same change trend. It indicates a slower heat absorption due to the smaller temperature gradient between the inlet and the artery blood temperature.

4.3. Particle deposition Fig. 13 plots the deposition fraction of 3μm heated particles, which is the same temperature as the inlet air, in various regions of the human upper airway model, with different inhalation flow rates under an inlet temperature of 80 oC. Clearly, the higher the inhalation flow rate, the more particles are deposited in the current model. The deposition fraction in the oral cavity is relatively high, taking up about one third of the total deposition. As in the oral cavity of the current model, the direction of the inlet hot air has a big angle with the figure inside, causing a large deposition on the part changing direction. The deposition in all regions increases with the breath intensity, with the most rapid growth in the pharynx-larynx region, which might be attributed to the turbulence flow field in the oropharynx. The deposition ratio in the trachea among the total deposition declines compared to that in the pharynx-larynx with the increasing flow rates. This is due to the better mixing and more uniform airstream in the trachea. As for the health implications, a low inhalation flow rate is suggested when faced with a fire, and the most likely region of thermal injury would be in the mucosa of the pharynx and larynx. The effects of the inhalation temperature are analyzed under the flow rate of 10L/min, seen in Fig. 14. With the inhalation the air temperature rises, the deposition fraction in the oral cavity, nasal-pharynx-larynx and the trachea are enhanced. The increasing deposition fraction in the human upper airway model is attributed to the thermophoretic force and the growing temperature gradient near the wall. Among the total deposition, the deposition proportion in the oral cavity increases from 1/3 to about 1/2 as the inlet temperature rises. It seems the heated particles with higher temperatures are more likely to deposit in the oral cavity, while the impact of the inhalation temperature gradually declines with the rising particle diameter. 23

Fig. 13. Deposition fraction of 3μm particles with different inhalation flow rates under the inlet air temperature of 80 oC.

Fig. 14. Deposition fraction of 3μm particles with different inhalation temperature under the inhalation flow rate of 10L/min. 24

Fig. 15 illustrates the effects of the particles’ diameters on the deposition rate. The deposition fraction in the oral cavity greatly increases with the particles’ diameter. However, the pharynx-larynx and trachea region witness a rising deposition from 3μm to 5μm and a reduction from 5μm to 6.5μm. The larger diameter particles are more likely to be deposited in the entrance and the lighter and smaller ones usually need longer lengths for deposition in the oral and pharynx-larynx region.

Fig. 15. Deposition fraction with different particles diameters under inhalation flow rate of 30L/min and the inlet air temperature of 80 oC.

5. Conclusion A realistic human upper airway model was created and experimentally validated, to simulate the steady oral inhalation of particles. It indicates that the oral cavity is more likely to be affected by the outside hot air. The mucosa in the oral, pharynx and larynx are more likely to incur thermal injury in the fire smoke situation. The inhalation flow rate is one of the main factors to affect the distribution of the interior temperature and particle deposition. With a higher inspiration flow rate, 25

the heated particles are more likely to be deposited, especially in the pharynx-larynx region. The increase of the inhalation air temperature and particles’ diameter give rise to a larger deposition fraction in the whole model. The heated particles with a higher temperature are more likely to deposit in the oral cavity, and the influence of the inlet temperature has only a minor influence in the case of a bigger particle diameter. In the further study, the effects of the human breathing pattern and mode (such as nasal and oral in some ratio of circulation breathing), particle aggregate and hygroscopicity, the heat transfer between the flow and the particles, and the evaporation in the mucous layer will be considered.

Acknowledgements We thank Prof. Jiyuan Tu and Dr. Kiao Inthavong for proving damping function UDF. This study was supported by the Chinese National Key Program for Research and Development (Grant No.2016YFC0802801 and 2016YFC0802807), National Natural Science Foundation of China (Grant No. 71473147), and the Open Foundation of Hefei Institute for Public Safety Research, Tsinghua University. The authors are deeply grateful to these supports.

References Ball, C. G., Uddin, M., Pollard, A., 2008. High resolution turbulence modelling of airflow in an idealized human extra-thoracic airway. Comput. Fluids. 37, 943-964. Cheng, Y. S., Zhou, Y., Chen, B. T., 1999. Particle deposition in a cast of human oral airways. Aerosol Sci. Tech. 31, 286-300. Desai, M. H., Rutan, R. L., Herndon, D. N., 1989. Managing Smoke Inhalation Injuries. Postgrad Med. 86, 69-76. Elcner, J., Lizal, F., Jedelsky, J., Jicha, M., Chovancova, M., 2016. Numerical investigation of inspiratory airflow in a realistic model of the human tracheobronchial airways and a comparison with experimental results. Biomech. Model Mechan. 15, 447-469. Goodarzi-Ardakani, V., Taeibi-Rahnia, M., Salimia, M.R., Ahmadi, G., 2016. Computational simulation of temperature and velocity distribution inhuman upper respiratory airway during inhalation of hot air. Respir. Physiol. Neurobiol. 223,49-58. Gosman, A. D., Loannides, E., 1983. Aspects of Computer Simulation of Liquid-Fueled Combustors. Journal of Energy. 7, 482-490. Grgic, B., Finlay, W. H., Burnell, P. K. P., Heenan, A. F., 2004. In vitro intersubject and 26

intrasubject deposition measurements in realistic mouth–throat geometries. J. Aerosol Sci. 35, 1025-1040. Grgic, B., Finlay, W. H., Heenan, A. F., 2004. Regional aerosol deposition and flow measurements in an idealized mouth and throat. J. Aerosol Sci. 35, 21-32. Heenan, A. F., Matida, E., Pollard, A., Finlay, W. H., 2003. Experimental measurements and computational modeling of the flow field in an idealized human oropharynx. Exp. Fluids. 35, 70-84. Hsu, D. J., Swift, D. L., 1999. The Measurements of Human Inhalability of Ultralarge Aerosols in Calm Air using Mannikins. J. Aerosol Sci. 30, 1331-1343 Inthavong, K., Ge, Q. J., Li, X. D., Tu, J. Y., 2012. Detailed predictions of particle aspiration affected by respiratory inhalation and airflow. Atmos. Environ. 62, 107-117. Inthavong, K., Tu, J.Y., Ahmadi, G., 2009a. Computational modelling of gas-particle flows with different particle morphology in the human nasal cavity. J. Comput. Multiphase Flows. 1, 57-82. Inthavong, K., Wen, J., Tian, Z. F., Tu, J. Y., 2008. Numerical study of fibre deposition in a human nasal cavity. J. Aerosol Sci. 39, 253-265. Inthavong, K., Wen, J., Tu, J.Y., Tian, Z.F., 2009b. From CT scans to CFD modelling fluid and heat transfer in a realistic human nasal cavity. Eng. Appl. Comp. Fluid. 3, 321-335. Jayaraju, S. T., Brouns, M., Lacor, C., Belkassem, B., Verbanck, S., 2008. Large eddy and detached eddy simulations of fluid flow and particle deposition in a human mouth–throat. J. Aerosol Sci. 39, 862-875. Jayaraju, S. T., Paiva, M., Brouns, M., Lacor, C., Verbanck, S., 2008. Contribution of upper airway geometry to convective mixing. J. Appl. Physiol. 105, 1733-1740. Johnstone, A., Uddin, M., Pollard, A., Heenan, A., Finlay, W. H., 2004. The flow inside an idealized form of the human extra-thoracic airway. Exp. Fluids. 37, 673-689. Kleinstreuer, C., Zhang, Z., 2010. Airflow and Particle Transport in the Human Respiratory System. Annu. Rev. Fluid Mech. 42, 301-334 Kobayashi, T., Sandberg, M., Kotani, H., Claesson, L., 2010. Experimental investigation and CFD analysis of cross-ventilated flow through single room detached house model. Build Environ. 45, 2723-2734. Lin, T. C., Breysse, P. N., Laube, B. L., Swift, D. L., 2001. Mouthpiece diameter affects deposition efficiency in cast models of the human oral airways. J. Aerosol Med. 14, 335-341. Lizal, F., Belka, M., Adam, J., Jedelsky, J., Jicha, M., 2015. A method for in vitro regional aerosol deposition measurement in a model of the human tracheobronchial tree by the positron emission tomography. Proc. IMechE, Part H: J. Engineering in Medicine, 229, 750-757. Lv, Y., Liu, J., Zhang, J., 2006. Theoretical evaluation of burns to the human respiratory tract due to inhalation of hot gas in the early stage of fires. Burns. 32, 436-446. Matida, E.A., Finlay, W.H., Lange, C.F., Grgic, B., 2004. Improved numerical simulation of aerosol deposition in an idealized mouth–throat. J. Aerosol Sci. 35,1-19. 27

Menter, F. R., Langtry, R., Völker, S., 2006. Transition Modelling for General Purpose CFD Codes. Flow, Turbulence and Combustion. 77, 277-303. Morsi, S. A., Alexander, A. J., 1972. An investigation of particle trajectories in two-phase flow systems. J. Fluid Mech. 55, 193. Park, S. S., Wexler, A. S., 2008. Size-dependent deposition of particles in the human lung at steady-state breathing. J. Aerosol Sci. 39, 266-276. Phuong, N. L., Ito, K., 2015. Investigation of flow pattern in upper human airway including oral and nasal inhalation by PIV and CFD. Build Environ. 94, 504-515. Phuong, N. L., Yamashita, M., Yoo, S.J., Ito, K., 2016. Prediction of convective heat transfer coefficient of human upper and lower airway surfaces in steady and unsteady breathing conditions. Build Environ. 100, 172-185. Rong, Y., Liu, W., Wang, C., Ning, F., Zhang, G., 2011. Temperature distribution in the upper airway after inhalation injury. Burns. 37, 1187-1191. Sandsund, M., Reinertsen, R. E., Holand, B., Bjermer, L., 2007. Thermoregulatory and respiratory responses in asthmatic and nonasthmatic subjects breathing cold and warm air during exercise in the cold. J. Therm Biol. 32, 246-254. Segal, R. A., Martonen, T. B., Kim, C. S., 2000. Comparison of computer simulations of total lung deposition to human subject data in healthy test subjects. J. Air Waste Manage. 50, 1262-1268. Stapleton, K. W., Guentsch, E., Hoskinson, M. K., Finlay, W. H., 2000. On the suitability of k-epsilon turbulence modeling for aerosol deposition in the mouth and throat: A comparison with experiment. J. Aerosol Sci. 31, 739-749. Tian, Z. F., Inthavong, K., Tu, J. Y., 2007. Deposition of inhaled wood dust in the nasal cavity. Inhal. Toxicol. 19, 1155-65. Wan, J., Zhang, G., Qiu, Y., Wen, C., Fu, T., 2016. Heat dissipation by blood circulation and airway tissue heat absorption in a canine model of inhalational thermal injury. Burns. 42, 548-555. Wang, S.M., Inthavong, K., Wen, J., Tu, J.Y., Xue, C.L., 2009. Comparison of micron- and nanoparticle deposition patterns in a realistic human nasal cavity. Resp. Physiol. Neurobi. 166, 142-151 Watt, P. W., Willmott, A. G. B., Maxwell, N. S., Smeeton, N. J., Watt, E., Richardson, A. J., 2016. Physiological and psychological responses in Fire Instructors to heat exposures. J. Therm. Biol. 58, 106-114. White, M. D., Cabanac, M., 1995. Respiratory heat-loss and core temperatures during submaximal exercise. J. Therm. Biol. 20, 489-496. Xu, C., Nielsen, P. V., Gong, G., Jensen, R. L., Liu, L., 2015. Influence of air stability and metabolic rate on exhaled flow. Indoor Air. 25, 198-209. You, K., Yang, H., Kym, D., Yoon, J., Haejun, Y., Cho, Y., Hur, J., Chun, W., Kim, J., 2014. Inhalation injury in burn patients: Establishing the link between diagnosis and prognosis. Burns. 40, 1470-1475. 28

Zhang, Z., Kleinstreuer, C., Kim, C. S., 2009. Comparison of analytical and CFD models with regard to micron particle deposition in a human 16-generation tracheobronchial airway model. J. Aerosol Sci. 40, 16-28. Zhang, Z., Kleinstreuer, C., 2003a. Species heat and mass transfer in a human upper airway model. Int J. Heat Mass Tran. 46, 4755-4768. Zhang, Z., Kleinstreuer, C., 2003b. Low-Reynolds-number turbulent flows in locally constricted conduits: A comparison study. Aiaa J. 41, 831-840. Zhang, Z., Kleinstreuer, C., Donohuec, J.F., Kim C.S., 2005. Comparison of micro- and nano-size particle depositions in a human upper airway model. J. Aerosol Sci. 36, 211-233. Zhang, Z., Kleinstreuer, C., 2011a. Laminar-to-turbulent fluid-nanoparticle dynamics simulations: Model comparisons and nanoparticle-deposition applications. Int. J. Numer. Meth. Bio. 27, 1930-1950. Zhang, Z., Kleinstreuer, C., 2011b. Computational analysis of airflow and nanoparticle deposition in a combined nasal–oral–tracheobronchial airway model. J. Aerosol Sci. 42, 174-194. Zhao, R., Di, L., Zhao, X., Wang, C., Zhang, G., 2013a. Measuring surface temperature and grading pathological changes of airway tissue in a canine model of inhalational thermal injury. Burns. 39, 767-775. Zhao, R., Di, L. N., Wen, C. Q., Ning, F. G., Zhang, G. A., 2013b. Circulational heat dissipation of upper airway: Canine model of inhalational thermal injury. Burns. 39, 1212-1220.

Highlights 

Simulation was conducted to investigate thermal injury in a human

airway model.



A three-dimensional realistic human upper airway model was built.



The SST transition turbulence model and a Lagrangian frame were

applied. 

Effects of flow rate, inhalation temperature and particle diameter

were studied. 

Oral cavity was proved to be more likely affected by the inlet air

conditions.

29