In this paper, the flow and heat transfer patterns in a subject-specific geometry of the human upper airway is numerically studied. The study was conducted for steady, inspiratory flow associated with quiet normal breathing with a tidal volume of V_{T} = 0.5 L/min and flow rate of Q = 250 cm^{3}/s. The numerical results confirmed in vivo measurement that the majority of heat transfer process takes place inside the nasal cavity. It is apparent that even for extreme cases (T_{∞} = −30 °C and T_{wall} = 37 °C), the inspired air approached the body temperature by the time it passes the distal nasopharyngeal region. The air temperature reached the body temperature by the time it is in the vicinity of the larynx.

## Introduction

The human lungs operate under a very tight physiological condition. Any small variation above or below its normal working condition makes the lungs susceptible to pathological conditions that impair normalcy. In order to keep the lungs healthy, inspired ambient air must be filtered and conditioned to the alveolar condition (i.e., fully saturated with water vapor at body temperature) before entering the lungs. In addition, potentially hazardous particles in the air must be prevented from entering the lungs. Studies by Webb [(1)], Cole [(2,3)], and Inglestedt [(4)] show that the heat and mass transfer associated with air conditioning process is primarily conducted in the nasal cavity. Many pathological conditions and interventional treatments of the upper respiratory system induce local geometrical changes that may affect the efficiency of the transport processes.

In the past, to understand the relation between airflow pattern and air conditioning in the nose and the effects of different pathological condition or surgical procedure, such as septal deviation, perforation or turbinates resection, experimental studies were conducted using in vivo measurement of temperature distribution in the human upper airways. However, in vivo measurement is restricted by two issues. The first issue is poor spatial and time resolution. The second one is the complexity of nasal cavity geometry severely limits access to entire nose. Advances in computer technology and numerical methods give rise to a number of numerical studies of the heat transfer process in human respiratory system. In general, these numerical studies can be grouped into two different classes. The first and most popular one considers that the conditioning of the inspired ambient air was carried out primarily in the nasal cavity. Therefore, most authors considered only nasal geometry in their study on respiratory transport phenomena [(5,6,7,8,9,10,11,12,13,14)]. The early published works on nasal cavity heat transfer employed simplified geometrical model and two-dimensional simulations at various transverse cross sections of the model [(12)]. The use of anatomically realistic models has been restricted by the ability of an imaging technique to capture the narrow passages in the nasal cavity and the meshing technique to generate mesh inside an extremely complex geometry. However, advances in medical imaging technology and the wide availability of commercial image processing and computational fluid dynamics (CFD) softwares alleviate some of the obstacles to use subject-specific nasal geometries. As a result, a number of works that employed subject-specific geometries, both healthy [(9,10,11,14)] and ones suffer from pathological conditions [(6,7,8),13], were published.

Another class of problems arises from the observation that there is a close relationship between respiratory heat and mass transfer and bronchoconstriction in patients suffering from exercise-induced asthma and the increase in risk of infection in patients with intubated airways. Only a limited number of publications belong to this class exist. They mainly concern with the intrathorachic region of the respiratory system and disregard the nasal cavity from their model. The most notable are the work of Ingenito et al. [(15)], Zhang and Kleinstreuer [(16)], and Tawhai and Hunter [(17)]. Ingenito et al., used finite difference methods to solve the heat transfer between the airway and the submucosal layer. Their airway model was based on a symmetrical dichotomously branching system. Zhang and Kleinstreuer [(16)] conducted a numerical simulation of heat and mass transfer in a human upper airway. The model of airways used in their study was a simplified human cast replica, consisting of a representative of extra-thoracic airway from oral cavity to trachea and a symmetric four-generation upper bronchial tree model. The steady 3D airflow is simulated for laminar as well as locally turbulent flow conditions using low Reynolds number (LRN) *k*–ω model. Using both symmetric and anatomically consistent airway model geometries Tawhai and Hunter simulated the thermofluid dynamics in the intubated airway from the trachea (Horsfield order 0) down to Horsfield order 6.

Most of published articles on respiratory heat transfer, in particular those that only consider nasal geometry, assume laminar flow condition due to low peak velocities and low Reynolds number [(11),12], though some turbulent conditions have been attempted [(6),13]. This assumption is supported by several experimental studies on nasal cavity flow which show that for normal breathing condition, airflow rate less than 180 cm^{3}/s on one side of cavity, the flow is predominantly laminar over most of the cavity even though some instability was experienced, they did not develop into a fully turbulent flow [(18),19]. However, if the domain of interest is beyond the nasopharyngeal region of the airway, even under normal breathing condition, the airflow distal to this region is highly turbulent with turbulent intensity ranges between 5% and 20%, as pointed by Lin et al. [(20)].

Thus, in this paper, the heat transfer in an anatomically realistic human upper airway geometry with the attached trachea and bifurcation is studied. Including the trachea is important to obtain realistic data both from flow and heat transfer perspective. The fluid flow model used here also includes a simple turbulence model. The remainder of the paper is organized as follows. Section 2 briefly describes the steps in the geometry reconstruction and the mesh generation processes. Section 3 presents the equations governing the flow and heat transfer and the numerical method used in this study. In Sec. 4, the results for both flow and heat transfer are presented. Finally, some conclusions are drawn in Sec. 5.

## Geometry Reconstruction and Mesh Generation

Historically, numerical studies on heat transfer processes in respiratory systems were carried out using simplified model of the airway geometries [(11),12,15,16,17]. Recent advances in medical imaging technology, digital image processing, and mesh generation techniques permit the numerical simulation of transport processes to be carried out on anatomically realistic geometries. These geometries were mostly reconstructed from a series of computed tomography (CT) images. However, high quality CT images are still an issue due to radiation dosage associated with increase in image resolution. As a result, the reconstruction of the airway geometry is the most laborious and time consuming step. In this study, the employed human upper airway geometry was reconstructed from pre-existing CT images of a patient with healthy airways. It comprises 390 slices of 512 × 512 pixels with 0.877 × 0.877 mm a spacing of 1.0 mm. The resolution of the scan is quite low. This is the normal standard for diagnosis in the United Kingdom. Owing to low resolution of the images, fully automatic methods such as 3D deformable model [(21)] can not be used in the segmentation process to delineate the airway boundaries.

The segmentation of the upper airway was carried out semimanually using a commercial software (amira). amira outputs a binary 3D array representing the pixels inside the reconstructed geometry. The surface mesh generation is the most complicated part of the preprocessing. We employed a standard marching cube type of algorithm to the tracheobronchial region and an advanced version of the marching cube method described in Ref. [22] was applied to the nasal cavity and nasopharynx. Before applying the advanced version of the marching cube method, the resolution of the nasal part of the geometry was doubled via linear interpolation. This resulted in better gradient approximation and control.

The unstructured volume tetrahedron mesh is generated using Delaunay triangulation method described in Refs. [23,24] prior to the application of different mesh improvement techniques described in Ref. [25]. In-depth discussion on techniques that was used to obtaine high quality mesh from a segmented geometry can be found in Ref. [26]. The final computational mesh for human upper airway geometry is composed of more than 2.1 × 10^{6} tetrahedra and the mesh is shown in Fig. 1.

## Governing Equations and Numerical Scheme

*ρ*is the density, $T\xaf$ is temperature,

*c*is the specific heat, and

_{p}*k*is the thermal conductivity of air. The mean laminar shear stress tensor is given as

_{T}is the turbulent eddy viscosity and δ

_{ij}is the Kronecker delta. In Spalart–Allmaras (SA) model, the eddy viscosity is a function of turbulent viscosity $\nu \u0302$ given by $\nu T=fv1\nu \u0302$ in which

*f*

_{v}_{2}is defined, respectively, as

*S*is the magnitude of vorticity. The parameter

*f*and the constant

_{w}*c*

_{w}_{1}are, respectively, given by

*c*

_{b}_{1}= 0.1355,

*σ*= 2/3,

*c*

_{b}_{2}= 0.622,

*κ*= 0.41,

*c*

_{w}_{2}= 0.3,

*c*

_{w}_{3}= 2 and

*c*

_{v}_{1}= 7.1.

By assuming that the physical properties of air are independent of temperature, the energy conservation equation can be decoupled from Eqs. 1,2 and uncoupled solution strategy can be adopted. This solution strategy involves obtaining the flow field first and followed by solving the energy conservation Eq. 3.

The solution for the flow problem was obtained using a class of finite element algorithm known as characteristic-based split (CBS) method [(27,28,29,30,31,32,33)]. The CBS scheme is very similar to the original Chorin split [(34),35] and also has similarities with the projection scheme widely employed in incompressible flow calculations. In order to solve the incompressible Navier–Stokes equations efficiently, the artificial compressibility (AC) method is used in conjunction with the CBS finite element method (FEM). Using the uncoupled solution strategy Eq. 3, in general, can be solved using any numerical scheme for solving convection-diffusion problem. In this study, we use explicit characteristic Galerkin FEM [(36)] for obtaining the temperature fields.

*β*is an artificial compressibility parameter.

The scheme then solves Eqs. 9,2 in three steps. In the first step, an intermediate velocity field is obtained, followed by the second step where pressure field is computed and, finally, the velocity field is corrected in the third step. The one equation SA turbulence model is added as a fourth step. The steps of the AC-CBS scheme in its semidiscrete form can be summarized as follows, omitting the overbar to simplify the notation

*Step 1: Intermediate momentum*

*Step 2: Pressure*

*Step 3: Momentum correction*

*Step 4: Transport of turbulence variable*

## Results and Discussion

The flow is assumed to be steady inspiratory flow for quiet normal breathing (tidal volume *V _{T}* = 0.5 L and flowrate

*Q*= 15 L/min). Since all the heat transfer aspects can be studied without loss of accuracy using steady state assumption, transient algorithm was not employed. Standard air properties are used in the calculations. Two cases are considered for this study. The first is associated with the condition where the inspired temperature is lower than the temperature of wall of the airways and the second one is where the inspired temperature is higher than the temperature of the wall of the airway. The temperatures of the inspired air are

*T*= −30 °C and

_{∞}*T*= 45 °C for the first and the second cases respectively. The temperature of the wall is assumed to be

_{∞}*T*

_{wall}= 37 °C and constant throughout the airway wall.

*x*

_{1}and

*x*

_{2}axes are in the plane of the inlet and the

*x*

_{3}co-ordinate is orthogonal to this plane then the equation for the velocity component normal to the plane reads as

*u*is the magnitude of velocity normal to the inlet plane and $\u2202p/\u2202x3|\Gamma inlet$ is the pressure gradient assumed to be constant across the inlet plane,

_{n}*μ*is the dynamic viscosity. The integration is performed by the FE method using the existing triangulation of inlet/outlets. In case of circular pipe, this method gives the parabolic Poiseuille velocity profile, but gives a natural generalization of that profile for more complicated inlet shapes. This way of implementing the inlet/exit boundary conditions clearly reduces the computation required for extension tubes. The outlet boundary condition is imposed by applying $p=0|\Gamma outlet$.

Figure 2 shows the stream traces of massless particles released from the right and left cavity, respectively. The flow is characterized by the formation of flow circulations at the anterior part of nasal vestibule and atrium of both left and right cavities and another one at the nasopharyngeal region. Other prominent flow features are that only very small proportion of air flow through olfactory slit and meatuses and majority of the flow is close to the nasal septum. This result is in agreement with the experimental results conducted by Kelly et al. [(37)] and Doorly and co-workers [(19)]. The only difference is that the circulation in the nasopharyngeal region did not appear in the experimental result. This is due to the fact the anatomy of every subject differs. Both experiments mentioned were conducted using a unilateral model of nasal cavity, while both cavities were included in the numerical model. The fact that there is always a certain degree of asymmetry between left and right cavity causing the flow to rotate and form a circulation in the area where the flow from both cavities meets.

To analyze the results, the velocity and temperature fields will be plotted on several coronal and transversal planes. The locations of these planes are depicted in Fig. 3. Plane C1, C2, C3, C4, C5, C6, C7, and C8 are located, respectively, at 1, 2, 3, 4, 5, 6, 7, and 8 cm from the most left point of the nasal geometry. For the transversal planes, planes T1, T2, T3, T4, T5, and T6 are located, respectively, at 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 cm from the lowest point at the nasal floor.

Temperature distribution plots for *T _{∞}* = −30 °C and

*T*= 45 °C are presented in Figs. 4,5. These figures demonstrate that the effectiveness of the nasal interior structures in heat exchange process. At steady state condition, the temperature of inferior region of the nasal cavity (T1 and T2 planes) is very close to wall temperature with exception of the nasal vestibule region, main nasal cavity and the distal nasopharyngeal region. In the middle region (T3, T4, and T5 planes), the area of low temperature in the left cavity extend further to the distal part of the nasopharyngeal region while in the right cavity the air temperature is very close to the wall temperature. Superior to this region (T6 plane) only the air in the nasal cavity and left superior meatus has temperature lower than the wall temperature. The air inside the sinuses has the same temperature as the wall.

_{∞}The minimum temperature distribution across different lines perpendicular to the coronal axis on the transversal planes for *T _{∞}* = −30 °C is shown in Fig. 6a. The graphs shows that on both cavity, the temperature of the inspired air rises sharply at a distance between 2 and 4 cm from the most anterior part of the nares that lays between the nasal valve and the anterior part of the middle meatus. This is the region in the nasal cavity where air velocity is high and the convection of heat is high. At the distal part of the nasopharyngeal region, the graph dip slightly as the result of air with low temperature, which flow through the middle and superior meatuses mix with the high temperature air that flow through both left and right inferior meatuses.

The Nusselt number distribution for cases where the ambient air temperature is *T _{∞}* = 45 °C and

*T*= −30 °C are shown in Figs. 7a and 7b, respectively. Both cases show a similar distribution pattern. The main heat transfer process happens in the nasal vestibules, the atrium of the right cavity, middle meatuses and the region near the septum. The rest of the upper airways is characterized by low Nusselt number where conduction dominates the heat transfer process. Higher air velocity and stronger recirculation flow inside the right cavity are responsible for the wider area of higher Nusselt number on the right cavity compare to the left one. However, the area of higher Nusselt number on the middle meatus of the left cavity is stretched further back.

_{∞}The assumption that the nasal cavity wall can supply/absorb all the required heat to bring the temperature of the inspired air to a near alveolar temperature is somewhat unrealistic. However, this issue is very complex. It will require the knowledge of the blood flow in the vascular bed in the nasal cavity wall and the heat transfer process from the epithelial tissue to the vascular bed. This knowledge is not available at the moment.

## Conclusions

The present study investigates the steady state heat transfer in subject-specific human upper airways. The subject-specific geometry of the upper airways was reconstructed from a series of CT images of resolutions normally employed for standard clinical diagnostic. We adopted the uncoupled solution strategy where the flow field is obtained first using the AC-CBS scheme followed by solving the thermal energy conservation equation using a fully explicit characteristic Galerkin FEM to obtained the temperature field.

The flow solution is able to capture qualitatively the main features of the flow through a healthy airway. The results of heat transfer calculation show that by the time inspired air reaches the nasopharyngeal region its temperature is within 10%–15% off the lung operating temperature. Distal to the larynx, we found that the air temperature is equal to the wall temperature. In both cases studied, the region between 2 and 4 cm from the most anterior part of the nares is the region with high temperature gradient. This region lay between the nasal valve and the anterior part of the middle meatus. Any change in the geometry of the nasal interior structure in this region is very likely to severely affect the heat transfer capacity of the nose.

This finding confirms the in vivo measurement that, under quiet normal breathing, the majority of the heat transfer processes is happening in the nasal cavity and shows the effectiveness of the nasal internal structures to condition the inspired air temperature closer to the lung operating temperature. In order to further improve the understanding of the respiratory heat and mass transfer the water transport needs to be taken into consideration. In addition, the influence higher flow rate and transient on respiratory heat and mass transfer will be considered in the future works.

## Acknowledgment

This work is supported by EPSRC Grant No. D070554.

## Nomenclature

*c*=_{p}specific heat at constant pressure

*k*=thermal conductivity

- $p\xaf$ =
mean pressure

*S*=magnitude of vorticity

*T*_{wall}=wall temperature

*T*=_{∞}ambient temperature

- $u\xafi$ =
mean velocity components

*β*=artificial compressibility parameter

*μ*=dynamic viscosity

*ν*=kinematic viscosity

*νT*=turbulent eddy viscosity

- $\nu \u0302$ =
turbulent viscosity

*ρ*=density

- $\tau \xafij$ =
mean laminar stress tensor components

- $\tau ijR$ =
Reynolds stress tensor components