CFD analysis of the flow structure in a monkey upper airway validated by PIV experiments

Respiratory Physiology & Neurobiology 271 (2020) 103304

Contents lists available at ScienceDirect

Respiratory Physiology & Neurobiology journal homepage: www.elsevier.com/locate/resphysiol

CFD analysis of the flow structure in a monkey upper airway validated by PIV experiments

T

Nguyen Lu Phuonga,b, , Tran Van Quangb, Nguyen Dang Khoac, Ji-Woong Kimd, Kazuhide Itoa ⁎

a

Faculty of Engineering Sciences, Kyushu University, Japan Faculty of Environment, University of Natural Resources and Environment, Hochiminh City, Viet Nam c Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, Japan d Korea Institute of Civil Engineering and Building Technology, Republic of Korea b

ARTICLE INFO

ABSTRACT

Keywords: Monkey upper airway Inhalation toxicology Computational fluid dynamics Particle image velocimetry Flow structure

Inhalation exposure to airborne contaminants has adverse effects on humans; however, related research is typically conducted using in vivo/in vitro tests on animals. Extrapolating the test results is complicated by anatomical and physiological differences between animals and humans and a lack of understanding of the transport mechanism inside their respective respiratory tracts. This study determined the detailed air-flow structure in the upper airway of a monkey. A steady computational fluid dynamics simulation, which was validated by previous particle image velocimetry measurements, was adopted for flow rates of 4 L/min and 10 L/min to analyze the flow structure from the nasal/oral cavities to the trachea region in a monkey airway model. The low Reynolds number type k–ε model provided a reasonably accurate prediction of the airflow in a monkey upper airway. Furthermore, it was confirmed that large velocity gradients were generated in the nasal vestibule and larynx regions, as well as increased turbulent air kinetic energy and wall sheer stress.

1. Introduction The rise in global and accessible travel has introduced a new class of infection risk through contact with new and unknown emerging pathogens. The control of airborne infection risks is particularly difficult because of the continuous intake of indoor/outdoor air through the human respiratory tract during steady breathing. Therefore, there is an urgent need for the prediction, control, and assessment of inhalation health risks. Understanding the transport mechanism inside the respiratory tract of humans and animals is essential for achieving a more comprehensive approach to inhalation toxicology research. Early in vivo experiments to investigate human nasal airflow patterns involved observations of the distribution of aerosol powder or smoke (Goodale, 1896; Parker, 1901). However, testing on the human nose produced a constriction or dilatation reaction in the nasal mucosa of the human subjects; swelling of the mucosa then affected nasal airflow patterns (Hahn et al., 1993). Therefore, in vitro techniques were introduced as an alternative method to avoid the ethical and medical problems related to human subjects. Nasal airflow has been examined using various types of in vitro measurement method. For example, laser Doppler velocimetry (LDV) has

been applied to measure the velocity fields of aerosolized water droplets in air through a human nasal model (Girardin et al., 1983). Moreover, Hornung et al. (1987) studied the flow of radioactive gas (xenon 133) through a nasal cavity model generated from a nasal cast of a cadaver. In contrast to these traditional forms of flow visualization; i.e., smoke, dye injection, hot wire anemometry, and LDV, particle image velocimetry (PIV) is considered a promising emerging method for analyzing flow dynamics, in which the 2D or 3D flow features can be quantified non-intrusively and instantaneously (Chung and Kim, 2008). Methods involving PIV have been applied to measure the velocity fields in physiological models of the upper human airway (Doorly et al., 2008; Hopkins et al., 2000; Kim and Haw, 2004; Phuong and Ito, 2015). To the best of our knowledge, our previous study (Kim et al., 2018) was the first to employ the PIV method to measure fluid flows within a 3D monkey airway model. The resulting PIV measurements yielded an in vitro 2D velocity map, which could provide useful data for the development and validation of numerical simulations. Recently, in silico studies using the computational fluid dynamics (CFD) technique have gained increasing attention because of their flexibility and lack of ethical constraints. CFD simulations are a promising technique for contributing to our understanding of the transport

⁎ Corresponding author. Present address: Interdisciplinary Graduate School of Engineering Sciences, Kyushu University, 6-1 Kasuga-koen, Kasuga, Fukuoka, 8168580, Japan. E-mail address: [email protected] (N.L. Phuong).

https://doi.org/10.1016/j.resp.2019.103304 Received 22 June 2019; Received in revised form 29 August 2019; Accepted 20 September 2019 Available online 20 September 2019 1569-9048/ © 2019 Elsevier B.V. All rights reserved.

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

mechanism of pollutants/pathogens in the respiratory tract. Previous simulations have determined detailed airflow patterns, air-conditioning, uptake of gas-phase contaminants, and particle deposition using realistic respiratory tract models for test animals and humans (Dong et al., 2018; Jayaraju et al., 2008; Kimbell et al., 1997; Lin et al., 2007; Morgan et al., 1991; Subramaniam et al., 1998; Wen et al., 2008). However, although a number of studies have simulated the airflow in the upper airway of humans and test animals, there is a lack of similar research for monkey airways. Kepler et al. (1998) developed a computer simulation model of a rhesus monkey to predict inspiratory nasal airflow and inhaled gas uptake. Their simulation results provided clear support for a hypothesis that inhaled airflow patterns have a significant impact on the location of primate nasal lesions induced by inhaled reactive gases. Moreover, the airflow, heat, and vapor exchange was simulated in the nasal passage of humans, chimpanzees, and macaques using a CFD method (Nishimura et al., 2016). As for rat, monkeys, and humans, Corley et al. (2012) reported the airflow and acrolein regional uptake within the respiratory tracts from the nose to lung branches, which are sensitive to airway geometry, airflow rates, acrolein concentrations, air:tissue partition coefficients, tissue thickness, and the maximum metabolism rate. Furthermore, comparative numerical models of inhaled micron-sized particle deposition in humans and animals have been the subject of sophisticated inhalation toxicological studies (Kabilan et al., 2016; Phuong et al., 2019; Shang et al., 2016). These previous studies have revealed the air conditioning function, heterogeneity of the gas-phase contaminant adsorption flux, and major particle deposition sites for various particle diameters through CFD analysis. Although airflow patterns, air conditioning, and gas-phase/ particle-phase contaminant transport mechanisms have been comprehensively studied by numerical analyses, CFD simulations have not been sufficiently validated by experimental results and we still lack a quantitative understanding of flow patterns in monkey upper airways. The nasal and oral cavities of animals are characterized by complex geometries and a wide range of target Reynolds numbers (Re) from laminar to turbulent; hence, the prediction accuracy and uncertainty of CFD simulation must be carefully validated using experimental results. This study aims to develop a realistic in silico model of a monkey upper airway for CFD analysis. CFD methods are advantageous as they can provide detailed 3D flow information that is typically difficult to obtain through in vivo/in vitro experiments. Variations in airway flow patterns and turbulence characteristics such as turbulent kinetic energy, wall shear stress, and static pressure distributions are also presented and discussed in this study.

Fig. 1. Schematic of the computational domain of the monkey airway.

Fig. 2. The surface mesh and a cross-sectional slice in the nasal cavity of the model and close-up of the prismatic boundary layer of the slice. Table 1 Details of CFD numerical condition setup. Numerical and boundary condition for CFD Turbulence model Mesh design Algorithm Scheme Inflow Outflow

2. Methods

Wall treatment Others

Laminar model SST k - ω model Low Reynolds Number type k - ε (Abe Nagano Kondo) model 7.6 million (unstructured) Steady state with SIMPLE algorithm Convection term: second order upwind Pressure inlet, Turbulent intensity(%) 10 Uout = -1.15 m/s (2 L/min) Uout = -2.3 m/s (4 L/min)Uout = -3.45 m/s (6 L/min) Uout = -4.6 m/s (8 L/min) Uout = 5.7 m/s (10 L/min) Velocity : no slip Isothermal condition

2.1. Segmentation of 3D model and mesh generation in our previous study; in this study, we adopted the exact same 3D geometrical data as previously reported visualization experiments (Ito et al., 2017; Kim and Phuong, 2018)

The subject monkey, Macaca fascicularis, was a 6-month-old male with a weight of 1.2 kg. Medical CT images were taken and converted into a compatible file format using Mimics ® (Materialise n.v.) to generate and modify 3D surface model of the medical images. A surface model was generated from continuous 2D contour data by translating segmented, modified, and smoothed contour points into a data series. 3D geometry model was then constructed from the segmentation by interpolation. External tubes for the nose and mouth were added to the entry of the airway from the nostril area, where connecting tubes provided the initiation airflow. The monkey respiratory tract model possessed a length and inner surface area corresponding to approximately 1.05 × 10−1 m and 2.81 × 10-5 m2, respectively. The monkey airway was divided into upper and lower airway regions as shown in Fig. 1. The 3D model was further optimized and exported in STL file format, which was loaded into the ANSYS 14.5 ICEM grid design module (ANSYS Inc., Canonsburg, PA). A hybrid grid was generated with grid refinement near the wall surfaces (Fig. 2). The details for generating monkey upper airway geometry have already been reported

2.2. Governing equations and turbulence models The conservation equations of mass and momentum were solved using the finite volume method. The continuity equation (Eq. (1)) and Reynolds-Averaged Navier–Stokes (RANS) equations (Eq. (2)) were numerically solved for steady, isothermal, incompressible, and turbulent flow:

(u¯ i ) =0 xi (u¯ i ) u¯ + u¯ j i = t xj

(1)

1 p + (2 Sij xi xj

u iu j)

(2)

where u¯ i is the ensemble-averaged velocity, p is the mean pressure, ρ 2

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

Based on the assumption that the Reynolds stress is proportional to the mean velocity gradient (Boussinesq hypothesis), the Reynolds stress equation can be written as Eq. (3):

u iu

j

= 2 t Sij

2 k 3

ij

(3)

where k is the turbulent quantity, ij is the Kronecker delta, and the turbulent or eddy viscosity (νt) is defined by eddy viscosity modeling (EVM) using RANS model. In recent years, different RANS turbulence models have been evaluated for transient-turbulent flows in the literature (Holbrook and Worth Longest, 2013; Jayaraju et al., 2008; Ma and Lutchen, 2009; Phuong and Ito, 2015, 2013; Stapleton et al., 2000; Zhang and Kleinstreuer, 2011a, 2011b). Holbrook and Worth Longest (2013) evaluated the performance of two-equation turbulence models in a respiratory double bifurcation model; the air flows were computed using laminar, standard k–ε, and low Reynolds number (LRN) k–ω turbulent models. Both the laminar and LRN k–ω models were validated by a clear correlation between particle deposition results and in vitro experimental data. Zhang and Kleinstreuer (2011a,b) conducted a Large Eddy Simulation (LES) and three RANS turbulence models: the LRN k–ω model, standard k–ω model, and shear stress transport (SST) k–ω transition model. The performances of LES, LRN k–ω, and SST k–ω transition models are not significantly different when predicting laminar flows and the transition to turbulent flow. Phuong and Ito (2015) compared the ability of four RANS turbulence models: the LRN Abe–Nagano–Kondoh (ANK) k–ε model, LRN Launder–Sharma (LS) model, renormalization group (RNG) k–ε model, and SST k–ω model, to predict laminar, transitional, and turbulent flows in a realistic human respiratory tract. The turbulence model LRN ANK showed the most agreement with PIV data in comparison to the other models at the maximum velocity location. Therefore, based on previous evaluations, the laminar model and two RANS models; i.e., a LRN ANK k–ε model and SST k–ω model, were adopted to predict the airflow field (Abe et al., 1994; Abe and Kondoh, 1995; Menter, 1994). 2.3. Boundary and initial conditions The transport equations of these models were discretized with a segregated solver using the SIMPLE (Semi-Implicit Method for PressureLinked Equations) algorithm to couple pressure and velocity, as well as second order upwind schemes for the momentum, turbulent kinetic energy, and specific dissipation rate. Air flow rates, Qin, of 4 L/min and 10 L/min were selected to pass through the airway model. The Reynolds numbers for steady flow measurement are approximately 900 for the 4 L/min air flow rate, and 2300 for the 10 L/min air flow rate. From the fluid dynamic aspect, flow rates (Qin) of 4 and 10 L/min represent laminar and transition regimes, respectively. Furthermore, Kim et al. (2018) reported that the flow not yet developed at flow rates of

Fig. 3. Contour of normalized velocity distributions of PIV and CFD under oral inhalation condition (Qin = 4 L/min) in the upper and lower airway regions (Un = um/ Ut , um = (u2 + v 2 ) ). PIV data obtained from Kim et al. (2018).

denotes the air density, u is the fluctuating velocity component, ν is the kinematic viscosity, and Sij is the mean rate-of-strain tensor. The ensemble-averaged equations can be solved if the unknown Reynolds stresses ( u i u j ) in Eq. (2) are related to the mean flow quantities.

Fig. 4. Profiles of normalized scalar velocity under oral inhalation condition with Qin = 4 L/min. PIV data obtained from Kim et al. (2018). 3

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

Fig. 5. A comparison two-dimensional normalized velocity vectors (Un) in central sagittal plane of oral inhalation condition in the upper airway region (a, b) and the lower airway region (c, d) with experimental data of Kim et al. (2018).

4 L/min. And the flow rate of 10 L/min was found to be the transitional flow regime. The previously reported PIV experiments were conducted under isothermal conditions; thus, the CFD analysis applied the same conditions. Flow fields in the nostril region were extended far upstream so that the velocity exhibited fully developed flow profiles with no significant radial velocity component in order to agree with the previous PIV experimental setup. The no-slip boundary condition was applied at the walls. The convergence criteria for the airflow was considered when the residuals for each of the equations and continuity equation dropped to 10−5 of the value. The CFD boundary and

numerical conditions are shown in Table 1. By applying SST k-ω and low Reynolds k-ε models, the wall unit y+ values were maintained on the order of less than 1 at all wall boundaries. The center of the prism grid cells nearest the wall surface was considered a non-dimensional distance y+, which was calculated using Eq. (4).

y+ =

uf yd

where yd is the distance normal to the wall surface, 4

(4) is the kinematic

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

Fig. 6. A comparison two-dimensional normalized velocity vectors (Un) in central sagittal plane of nasal inhalation condition in the upper airway region (a, b) and the lower airway region (c, d) with experimental data of Kim et al. (2018).

viscosity, and uf is the friction velocity. The friction velocity is defined as

uf =

w

3. Results and discussion 3.1. Turbulent model selection

(5)

To evaluate the occurrence of laminar, transitional, and turbulent flows in the monkey upper airway model, the laminar, SST k-ω, and LRN ANK k-ε models are considered for a 4-L/min flow rate. We compare the flow patterns (qualitative calculation) and flow profiles (quantitative calculation) obtained by the three models. PIV data are utilized in this study obtained from previous study, Kim and Phuong,

where τw is the wall shear stress. The grid independence was also carefully investigated and the solutions obtained with grid sizes of 3,250,000, 5,450,000, 7,620,000, and 8,970,000 control volumes were also investigated (Appendix A). Finally, we adopted 7,620,000 elements for the CFD analysis in this study. 5

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

Fig. 7. A comparison profiles of normalized velocity under oral inhalation condition with PIV data of Kim et al. (2018).

2018. Normalized velocity distributions of PIV and CFD are selected at the same plane location and shown in Fig. 3. In order to compare the flow pattern of two-dimensional PIV experiments and CFD predictions, the 2D normalized scalar velocity (Un) was defined as the mean velocity (um) at the central sagittal plane divided by the velocity at the trachea (Ut) (Un = um/Ut, um = u2 + v 2 , where u = velocity magnitude of xcomponent, v = velocity magnitude of y-component). Fluid flows in the PIV and CFD results exhibited a similar flow structure in the upper region and complex flow patterns in the lower region. The angle of the inflow boundary leads to high flow velocity in the vicinity of the upper wall of the oral cavity. After the fluid flow passes the larynx, the fluid velocity increases, and flow separation occurs with abrupt geometrical changes from the larynx to the trachea. Thus, the laminar, SST k-ω, and LRN k-ε turbulence models effectively capture the flow pattern in the

entire model. Specifically, shearing flow after the glottis is observed when using the laminar, low Reynolds k-ε, and SST k-ω models. An asymmetric laryngeal jet forms in the center of the larynx region and is deflected towards the anterior trachea. Downstream of the glottis, the magnitude of the tracheal jet is well predicted for the two peaks. In all cases, the scalar velocities are relatively low in the upstream, whereas the scalar velocity increases with distance downstream and is highest in the tracheal jet from the larynx to the trachea. In general, all models exhibit qualitatively good predictions of the air flow profile in the monkey upper airway. However, it is still important to validate the velocity profiles through a quantitative assessment using PIV values. Fig. 4a and b show the scalar velocity profile (Un) plotted at the vertical cross lines in the upper and lower airway regions for an oral inhalation flow rate of 4 L/min, respectively. The length of the cross 6

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

Fig. 8. A comparison profiles of normalized velocity under nasal inhalation condition with PIV data of Kim et al. (2018).

line (y) was normalized by the tracheal diameter (D). From the Fig. 4a, the peak velocity is found in cross line 2 when y/D had a value of between 0.5 and 1 for both the PIV experiment and CFD simulation. The flows in both the PIV measurements and in the modeled case in the oral cavity are in good agreement. The difference ratio with the peak values for the PIV experiment are about 4% according to the LRN k-ε model, 6% according to the Laminar and the SST k-ω models. Moreover, the maximum velocity fluctuations occur after the larynx and at the entrance of the trachea. Therefore, cross line 5 at the beginning of the trachea is selected to present the velocity profile in the lower airway region (see supplementary drawing Fig. 4b). Clearly, the simulated velocity profiles are almost identical for different simulated models. The flow trends for three models are similar to the PIV results; however, the predicted velocity magnitudes of the peak flow fluctuations vary

between the turbulence models. Both SST k-ω and LRN k-ε turbulence models can predict the peaks, as shown in Fig. 4b. The laminar model does not exhibit a good quantitative prediction of the velocity profile in this case. The velocity profiles are similar in the anterior trachea (first peak) but differ in the central trachea (second peak). The first peak velocities are observed at y/D between 0 and 0.25; the discrepancies between the peak values of the simulation models and the experimental data is 17%, 14%, and 8% for the laminar model, SST k – ω model, and LRN k – ε model, respectively. All three models overpredict the second peak velocities at a y/D of 0.4–0.5. The results of the LRN k-ε model are more consistent with the PIV data than those of the Laminar and SST kω models due to distinct differences in the velocity profile within the entire airway. As a result of this preliminary study, the LRN k-ε turbulence model (ANK model) was selected for further simulations. 7

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

Fig. 9. Normalized velocity distribution in cross sections.

3.2. Velocity vector map at the central sagittal plane

anterior wall to the end of the trachea. Formation of the laryngeal jet is also observed in human airways (Heenan et al., 2003; Jayaraju et al., 2008). The Reynolds number based on the inlet ranges from approximately 900 to 2300 for flow rates of 4–10 L/min. In Fig. 5c and d, the flows are characterized by a region of separated flow at the larynx and down to the trachea for both CFD and PIV (lower airway region). The separation and acceleration of flow is caused by the inclination and contraction of the cross-sectional area of the epiglottis and larynx. The high-speed flow enters the trachea at an angle, impinging against and flowing along the anterior wall of the trachea.

3.2.1. Oral inhalation condition The CFD results divided by upper and lower airway regions are compared with the corresponding vector map measured by PIV at inhalation flow rates of 4 L/min and 10 L/min (Fig. 5). The flow pattern through the oral cavity was obtained by PIV (left figure) and CFD (right figure). After entering the mouth, the recirculation zone, back flow, and reattachment point are observed in Fig. 5a and b (upper airway region). Reversed flow is confirmed as a result of the adverse pressure gradient caused by the abrupt increase in the cross-sectional area from the mouth opening to the oral passage. The fluid flow follows the upper part of the oral cavity until it turns into the pipe elbow (pharynx region) in two cases. The change in cross-sectional area of the larynx causes the flow to accelerate, resulting in a jet. The laryngeal jet follows the inner

3.2.2. Nasal inhalation conditions Velocity vectors corresponding to the nasal inhalation condition at the central sagittal plane are depicted in Fig. 6. In Fig. 6a and b, the normalized velocity vectors of the left nasal cavity provide qualitative 8

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

scalar velocities (normalized velocity) are less than 1 in the oral cavity region (Fig. 7a and b) and greater than 1 in the lower region (Fig. 7c). In Fig. 7a, the distribution of the normalized stream wise velocity component (u/Un) shows partly negative velocities in the vicinity of inlet opening (Line 1). In Fig. 7b, the normalized scalar velocity profile (Un) of Line 2 indicates that the fluid on the upper (vicinal palate) side moves faster than that on the lower (vicinal tongue) side. Hence, the two peaks occur at the trachea (Fig. 7c, Qin = 4 L/min). The results in Fig. 7c show the fluid on the anterior side moving faster while the fluid on the posterior side is slower. CFD predictions match the PIV data well at both inhalation flow rates via the oral cavity. A clear difference can be seen between the profiles for the lower (4 L/min) and higher (10 L/ min) flow rates in Fig. 7c. At the higher flow rate, the CFD velocity peak occurs in the middle region while the peak according to PIV data is on the anterior side. As shown in Fig. 8a (left nasal cavity), a pronounced concave curve profile is plotted with the protruding portion along Line 3. Fig. 8b (right nasal cavity) shows unstable profiles as a result of complex shapes. The concave curves in Fig. 8a and b are caused by the unique shape of the monkey and the contracting cross-sectional area in the nasal cavity. Similar to the results in Fig. 7c for Qin = 4 L/min, the fluid flow rate of 4 L/min in Fig. 8c is unevenly developed for both PIV and CFD results. When the oral inhalation flow rate is increased to 10 L/min, the velocity profile becomes well developed (Fig. 8c). In general, the CFD predictions are reasonably consistent with the results of PIV experimental data. The numerical results presented here have some similarities/consistencies from previous PIV experiments. The effect of changing Reynolds number can be seen in the normalized scalar velocity profiles plots of Figs. 7 and 8. The normalized velocity contour plots are quite similar for two air flow rates in the upper airway region, indicating an independent of flow to Reynold effects. However, there are some structural variations with increasing Reynolds number in the lower airway region. In case of higher Reynolds number, the curvature profile (Figs. 7c and 8 c) becomes more pronounced. 3.4. Velocity contour in the 2D cross sections The scalar velocity distribution at different cross sections of the airway as airflow enters the nostrils during inhalation conditions are presented in Fig. 9. The bulk of the flow enters through the left and right nostril, as shown in Fig. 9a and b, corresponding to flow rates of 4 L/min and 10 L/min, respectively. The bulk flows move in the upward direction at cross section A. The airflow changes direction to near the septum wall at cross section B (nasal vestibule region) and is influenced by the narrowed region. The velocity magnitude of the left nasal cavity becomes dominant over the right nasal cavity at a high flow rate (Qin = 10 L/min, Fig. 9b). At cross section D (turbinate region), the main flow is still close to the septum wall and low flow velocity occurs in the outer meatus region (Fig. 9a and b). Air flow then moves into the region of the nasal cavity with maxillary sinuses (slice E) and a stagnant region occurs in the sinus region of the airway (Fig. 9a and b). The flow across the epiglottis (slice H) is particularly concentrated in the center of the airway. The velocities are higher on the anterior side of the trachea, as shown by the contour plot across the larynx section (slice I). The flow then impacts the front of the tracheal wall and creates a high velocity area as the larynx jet moves deeper down the trachea in slice K (Fig. 9a). In contrast, the slice K results exhibit different features due to the higher nasal inhalation flow rate (Fig. 9b). A near uniform velocity profile is observed as a fully developed turbulent flow through the circular trachea (Fig. 9b, slice L). The air flow distribution in the upper region (slices A–G) is quite similar to those in a human nasal cavity reported by Inthavong et al. (2009). Fig. 9c and d show the scalar velocity contours at different cross sections in the airway during oral inhalation. As the flow reaches the beginning of the oral cavity, it is concentrated in the medial part (slice

Fig. 10. Turbulence Kinetic Energy (TKE) magnitude contour in the central sagittal plane and in horizontal planes for nasal and oral inhalation conditions.

visualization of the flow field. The normalized velocity vectors in the left nasal cavity at a flow rate of 4 L/min show flow separation immediately posterior to the narrowest valve (Fig. 6a). As the flow rate increases to 10 L/min (Fig. 6b), the flow moves upward to the middle turbinate region. Note that the human airway structure has three turbinates (inferior, middle, and superior) whereas a monkey airway has only two turbinates corresponding to the human inferior and middle turbinates (Chamanza and Wright, 2015). The scalar velocities are relatively low in the upstream, but the scalar velocity increases downstream and is highest in the laryngeal jet from the larynx to trachea at both inhalation flow rates (Fig. 6c and d). 3.3. Velocity profile in representative cross sections A quantitative comparison between the experimental and simulated velocity profiles of oral and nasal inhalation are shown in Figs. 7 and 8, respectively. The cross lines shown in the graph are indicated in the supplemental drawing in the top margin of Fig. 7. In Fig. 7, the peak 9

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

Fig. 11. Surface wall shear stress distribution mapping at Qin = 4 L/min.

A1) at both flow rates. However, the velocities are higher on the upper side of the oral cavity (slice B1). When the bulk flow reaches the oropharynx and larynx (slice F1 and G1), the airflow velocity reduces due to widening of the oral cavity. When the flow reaches the epiglottis, the bulk of the flow moves to the center (slice I). Clearly, the velocity moves to the anterior side due to area expansion and the transition to turbulent flow with large-scale momentum transfer. For oral inhalation, the scalar velocities at slices J, K, and L are much larger than those for nasal inhalation.

In this study, we observed a region of high velocity when the air flow passes through the mouth opening and travels above the tongue. Hence, the tongue is recognized as the largest obstacle to flow and aerosol deposition in human airway models (Grgic et al., 2004; Heenan et al., 2003). The location of the tongue in a human airway might affect the turbulent flow, as reported by Lin et al. (2007). Therefore, the tongue location in the monkey airway might also change the flow pattern from slice A1 to G1 in Fig. 9c and d. 10

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

3.6. Wall-shear stress distribution on the airway wall The effect of airflow on inducing wall shear stress has not yet been quantitatively determined for the monkey airway wall. Understanding wall-shear stress is useful for comprehending the resistive forces, exchange processes, near-wall flow field, and potential sites of particle deposition during respiration in humans (Elad et al., 2006; Gambaruto et al., 2010; Tian et al., 2017). Fig. 11 shows the wall-shear stress distributions with plots of the surface integrated wall-shear stress under nasal and oral inhalation conditions at Qin = 4 L/min. In Fig. 11a, the highest wall shear stresses are found in the vestibule region. The wall shear stress decreases gradually as the airway widens in the middle region. The increasing wall shear stress in the nasopharynx is caused by joining of the two cavities and a change in the flow direction whereby air is transported towards the trachea. The highest wall-shear stress distributions are observed in the larynx and trachea regions in Fig. 11a and b. As expected, the largest flow resistance occurs immediately after the epiglottis. The larynx jet impinges on the anterior protrusion of the larynx, creating a local maximum of wall-shear stress visible under both nasal and oral conditions (see anterior view in Fig. 11a and b). Simulations of the nasal inhalation process reveal wall shear stresses as high as 0.3 Pa in the upper region and 1.5 Pa in the lower region (larynx section). Similarly, during oral inhalation, the wall shear stresses are also 0.3 Pa in the upper region and 1.8 Pa in the lower region. In human airways, maximum shear stresses have been found in the range of 0.2 Pa in the nasal cavity at an inspiration rate of 20 L/min (Elad et al., 2006) and 0.3 Pa at a flow rate of 10 L/min (Inthavong et al., 2014). Increasingly high shear stresses in the concentrated regions may also cause irritation of the blood vessels and nasal lesions with exposure to air pollutants (Inthavong et al., 2014; Kepler et al., 1998).

Fig. 12. Comparison of pressure drop across nasal cavity as a function of flow rate during nasal inhalation condition.

3.5. Turbulent kinetic energy The contours of turbulent kinetic energy (TKE) under nasal inhalation conditions and in three cross sections for two flow rates are shown in Fig. 10a and b. In both cases, the region of high TKE is located in the trachea near the end of the jet core and is attributable to transition of the larynx jet to turbulence. TKE is reduced with distance downstream to the lower part of the trachea. The maximum values recorded are 0.2 m2/s2 and 5.5 m2/s2 for inhaled flow rates of 4 L/min and 10 L/min, respectively. TKE occurs in the region from the nostrils to the nasal cavity at Qin = 10 L/min. According to Fig. 10b, cross-sectional slices of higher flow rate (10 L/min) exhibit broader TKE values than those at lower flow rate (4 L/min). Fig. 10c and d show the contours of TKE in oral inhalation mode. The TKE values are observed in the lower airway region and reach a maximum at the entrance of the trachea. The increase in TKE can thus be linked to the increase in fluid velocity following the larynx. The cross section at slice J shows that strong turbulence activity occurs in the anterior side of the trachea. In the trachea region, maximum TKE values of 0.4 m2/s2 and 5.5 m2/s2 are observed for inhaled flow rates of 4 L/ min and 10 L/min, respectively. Regions of maximum TKE are recorded in slice J and slice K, but dissipate gradually at slice L (Fig. 10a–d). The maximum TKE of oral inhalation is higher than the maximum TKE in the nasal simulation. Overall, TKE increases with increasing inhalation flow rates.

3.7. Pressure difference Numerical predictions of pressure loss (pressure drop) are compared to in vitro measurement data available from Kelly et al. (2005) in Fig. 12. Inhalation flow rates of 2, 6, and 8 L/min have been added to simulate the pressure loss (ΔP) in the monkey upper airway. For nasal inhalation, the maximum pressure loss values are 32 Pa and 120 Pa for Qin = 4 L/min and 10 L/min, respectively. The pressure loss at different flow rates shows a non-linear trend and is in agreement within the range 1–10 L/min. Fig. 13 shows a comparison of area-averaged static pressure profiles of different cross-sections for nasal and oral inhalation conditions at flow rates of 4 L/min and 10 L/min. Less resistance is found during oral inhalation where the static pressure decrease is smaller than that during nasal inhalation. In the lower region of approximately 8 cm from the tip of the nose/mouth, significant airway resistance to airflow is observed for all cases. This large static pressure loss occurs in the larynx region where the nasal and oral breathing flows merge.

Fig. 13. Comparison of area-averaged static pressure on different cross-sections at flow rates of 4 L/min and 10 L/min. 11

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

3.8. Implications for the fluid structure in a monkey upper airway

models are based on solving for time-averaged equations, the computational cost to perform a simulation is quite reasonable. The RANS models often involve simplified assumptions in the inflow boundary condition, reasonable grid design, and steady-state analysis. Hence, the RANS turbulence models are the most preferred method for turbulent flow simulations if the prediction of mean flow characteristic is of main interest. In this study, we focus on testing the validity of the commonly RANS modeling approach, and saving the LES for the description of fluid behavior in the monkey airway models in the next stage.

The results of this study indicate the importance of correct inlet boundary conditions and the need to consider upstream effects in experimental and computational studies of complex monkey upper airway models. This study represents one of the first steps towards validating CFD data of a monkey upper airway using PIV measurement results. We reveal that oral inhalation and nasal inhalation in the upper region result in different flow structures in the lower region of the upper airway. Thus, the effect of turbulence on the airflow in different regions of the monkey upper airway is significant. The presence of large fluctuations (turbulent characteristics) in the vestibule and constriction region (epiglottis) increased the localized wall shear stress. CFD analysis revealed that regions of local maximum TKE and wall shear stress are associated with turbulent jet flow in the larynx region. The pressure loss increased with increasing inhalation flow rates. Moreover, the results of this study imply that CFD simulations can be effectively applied to precisely predict flow patterns and profiles and provide further detailed information of fluid behavior such as turbulent kinetic energy, wall shear stress, pressure loss, and static pressure, which cannot be derived using the PIV method. However, the absence of breathing cycles and cyclic wall motion in the model is a major limitation in the simulation of breathing flows. Nevertheless, the rigid-walled model still provides insights into basic flow properties. Further studies should also perform simulations under transient flow conditions. The selection of appropriate CFD models is critical for the validation process of monkey airway airflow with PIV measurement results. LES is more accurate than the RANS models since the large eddies of turbulent flow are computed directly and only sub-grid scale (small scale) motions are modelled in LES approach (Jayaraju et al., 2008; Liu et al., 2007). However, LES has some limitations to become a widely practical tool. The computational cost of LES has been reported that a few hundred times higher than that of RANS models (Li et al., 2017; Zhang and Kleinstreuer, 2011a, 2011b), and LES has required users with extensive years of experience on LES techniques. Whereas RANS turbulent

4. Conclusions This study performed numerical validation of a realistic monkey airway model incorporating adequate geometrical characteristics of upper airway structures. Particle image velocimetry results were used for the validation and discrepancies were analyzed between the experimental and CFD simulation results. The air flow calculation was very challenging due to the complex geometry of the model, characterized by narrow and bending passages. However, our results suggest that flow can be reasonably predicted in a complex model by the CFD technique with the LRN type k- ε model. This validation work can be used for the development of in silico models to complement and strengthen in vivo and in vitro experiments. Computational simulations have great potential to contribute to studies of therapeutic drug delivery and inhalation toxicology for nonhuman primate targets. Therefore, this study represents an important step towards extrapolating respiratory mechanisms from surrogate animals to humans for therapeutic or toxicological analyses. Acknowledgement The authors gratefully acknowledge the financial support provided by a Grant-in-Aid for Scientific Research (JSPS KAKENHI No 17F17078 and No 18H03807).

Appendix A Grid independence The computational model was discretized with unstructured tetrahedral elements. Prism layers were applied in the vicinity of the wall surfaces to analyze the precise flow profile inside the viscous sub-layers. The adequacy of the grid resolution was tested by refining the surface elements of the airway wall and the outlet using element sizes of 0.7 mm, 0.6 mm, 0.5 mm, and 0.4 mm in the monkey airway model. The total numbers of 3D elements generated according to these surface element sizes were 3,250,000, 5,450,000, 7,620,000, and 8,970,000, respectively. Fig. A1 shows the

Fig. A1. Grid independence testing for oral inhalation based on velocity profile at Line 6. 12

Respiratory Physiology & Neurobiology 271 (2020) 103304

N.L. Phuong, et al.

CFD results of the grid independence check with scalar velocities along Line 6. It was confirmed that there was critical discrepancy between the results using the four types of grid designs. To compromise between the computational cost and accuracy of the simulations, a hybrid mesh size of 7.6 million was used to discretize the model.

Kelly, J.T., Asgharian, B., Wong, B.A., 2005. Inertial particle deposition in a monkey nasal mold compared with that in human nasal replicas. Inhal. Toxicol. 17, 823–830. https://doi.org/10.1080/08958370500241270. Kepler, G.M., Richardson, R.B., Morgan, K.T., Kimbell, J.S., 1998. Computer simulation of inspiratory nasal airflow and inhaled gas uptake in a rhesus monkey. Toxicol. Appl. Pharmacol. 150, 1–11. https://doi.org/10.1006/taap.1997.8350. Kim, J.W., Phuong, N.L., Aramaki, S., ichiro Ito, K., 2018. Flow visualization through particle image velocimetry in realistic model of rhesus monkey’s upper airway. Respir. Physiol. Neurobiol. 251, 16–27. https://doi.org/10.1016/j.resp.2018.02.007. Kim, S.K., Haw, J.R., 2004. An investigation on airflow in pathological nasal airway by PIV. J. Vis. 7, 341–348. https://doi.org/10.1007/BF03181538. Kimbell, J.S., Godo, M.N., Gross, E.A., Joyner, D.R., Richardson, R.B., Morgan, K.T., 1997. Computer simulation of inspiratory airflow in all regions of the F344 rat nasal passages. Toxicol. Appl. Pharmacol. 145, 388–398. https://doi.org/10.1006/taap.1997. 8206. Li, C., Jiang, J., Dong, H., Zhao, K., 2017. Computational modeling and validation of human nasal airflow under various breathing conditions. J. Biomech. 64, 59–68. https://doi.org/10.1016/j.jbiomech.2017.08.031. Lin, C.-L., Tawhai, M.H., McLennan, G., Hoffman, Ea, 2007. Characteristics of the turbulent laryngeal jet and its effect on airflow in the human intra-thoracic airways. Respir. Physiol. Neurobiol. 157, 295–309. https://doi.org/10.1016/j.resp.2007.02. 006. Liu, Y., Matida, Ea., Gu, J., Johnson, M.R., 2007. Numerical simulation of aerosol deposition in a 3-D human nasal cavity using RANS, RANS/EIM, and LES. J. Aerosol Sci. 38, 683–700. https://doi.org/10.1016/j.jaerosci.2007.05.003. Ma, B., Lutchen, K.R., 2009. CFD simulation of aerosol deposition in an anatomically based human large-medium airway model. Ann. Biomed. Eng. 37, 271–285. https:// doi.org/10.1007/s10439-008-9620-y. Menter, F.R., 1994. Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications. pp. 32. Morgan, K.T., Kimbell, J.S., Monticello, T.M., Patra, A.L., Fleishman, A., 1991. Studies of inspiratory airflow patterns in the nasal passages of the F344 rat and rhesus monkey using nasal molds: Relevance to formaldehyde toxicity. Toxicol. Appl. Pharmacol. 110, 223–240. https://doi.org/10.1016/S0041-008X(05)80005-5. Nishimura, T., Mori, F., Hanida, S., Kumahata, K., Ishikawa, S., Samarat, K., MiyabeNishiwaki, T., Hayashi, M., Tomonaga, M., Suzuki, J., Matsuzawa, Tetsuro, Matsuzawa, Teruo, 2016. Impaired air conditioning within the nasal cavity in flatfaced Homo. PLoS Comput. Biol. 12, 1–18. https://doi.org/10.1371/journal.pcbi. 1004807. Parker, C.A., 1901. Some observations and remarks on the air-currents in nasal respiration. J. Laryngol. Rhinol. Otol. 16, 345–355. https://doi.org/10.1017/ S1755146300169690. Phuong, N., Khoa, N., Inthavong, K., Ito, K., 2019. Particle and inhalation exposure in human and monkey computational airway models. Inhal. Toxicol. https://doi.org/ 10.1080/08958378.2018.1545810. 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. https://doi.org/10.1016/j.buildenv.2015.10.002. Phuong, N.L., Ito, K., 2013. Experimental and numerical study of airflow pattern and particle dispersion in a vertical ventilation duct. Build. Environ. 59, 466–481. https://doi.org/10.1016/j.buildenv.2012.09.014. Shang, Y., Dong, J., Inthavong, K., Tu, J., 2016. Comparative numerical modeling of inhaled nanoparticle deposition in human and rat nasal cavities. Toxicol. Sci. 152, 284–296. https://doi.org/10.1093/toxsci/kfw087. Stapleton, K., Guentsch, E., Hoskinson, M., Finlay, W., 2000. On the suitability of k–ε turbulence modeling for aerosol deposition in the mouth and throat: a comparison with experiment. J. Aerosol Sci. 31, 739–749. https://doi.org/10.1016/S00218502(99)00547-9. Subramaniam, R.P., Richardson, R.B., Morgan, K.T., Kimbell, J.S., Guilmette, R.A., 1998. Computational fluid dynamics simulations of inspiratory airflow in the human nose and nasopharynx. Inhal. Toxicol. 10, 91–120. https://doi.org/10.1080/ 089583798197772. Tian, L., Shang, Y., Dong, J., Inthavong, K., Tu, J., 2017. Human nasal olfactory deposition of inhaled nanoparticles at low to moderate breathing rate. J. Aerosol Sci. 113, 189–200. https://doi.org/10.1016/j.jaerosci.2017.08.006. Wen, J., Inthavong, K., Tu, J., Wang, S., 2008. Numerical simulations for detailed airflow dynamics in a human nasal cavity. Respir. Physiol. Neurobiol. 161, 125–135. https:// doi.org/10.1016/j.resp.2008.01.012. Zhang, Z., Kleinstreuer, C., 2011a. Computational analysis of airflow and nanoparticle deposition in a combined nasal–oral–tracheobronchial airway model. J. Aerosol Sci. 42, 174–194. https://doi.org/10.1016/j.jaerosci.2011.01.001. Zhang, Zhe, Kleinstreuer, Clement, 2011b. Laminar-to-turbulent fluid–nanoparticle dynamics simulations: Model comparisons and nanoparticle-deposition applications. Int. J. Numer. Method. Biomed. Eng. 27, 1930–1950. https://doi.org/10.1002/cnm.

References Abe, K., Kondoh, T., 1995. ˜ P E R G a M O N a New Turbulence Model for Predicting Fluid Flow and Heat Transfer in Separating and Reattaching Flows I1. Thermal Field Calculations. Abe, K., Kondoh, T., Nagano, Y., 1994. A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching Flow field calculations. Int. J. Heat Mass Transf. 37, 139–151. Chamanza, R., Wright, J.A., 2015. A Review of the Comparative Anatomy, Histology, Physiology and Pathology of the Nasal Cavity of Rats, Mice, Dogs and Non-human Primates. Relevance to Inhalation Toxicology and Human Health Risk Assessment. J. Comp. Pathol. 153, 287–314. https://doi.org/10.1016/j.jcpa.2015.08.009. Chung, S.-K., Kim, S.K., 2008. Digital particle image velocimetry studies of nasal airflow. Respir. Physiol. Neurobiol. 163, 111–120. https://doi.org/10.1016/j.resp.2008.07. 023. Corley, R.A., Kabilan, S., Kuprat, A.P., Carson, J.P., Minard, K.R., Jacob, R.E., Timchalk, C., Glenny, R., Pipavath, S., Cox, T., Wallis, C.D., Larson, R.F., Fanucchi, M.V., Postlethwait, E.M., Einstein, D.R., 2012. Comparative computational modeling of airflows and vapor dosimetry in the respiratory tracts of rat, monkey, and human. Toxicol. Sci. 128, 500–516. https://doi.org/10.1093/toxsci/kfs168. Dong, J., Shang, Y., Tian, L., Inthavong, K., Tu, J., 2018. Detailed deposition analysis of inertial and diffusive particles in a rat nasal passage. Inhal. Toxicol. 1–11. https://doi. org/10.1080/08958378.2018.1439549. 0. Doorly, D., Taylor, D.J., Franke, P., Schroter, R.C., 2008. Experimental investigation of nasal airflow. Proc. Inst. Mech. Eng. Part H J. Eng. Med. 222, 439–453. https://doi. org/10.1243/09544119JEIM330. Elad, D., Naftali, S., Rosenfeld, M., Wolf, M., 2006. Physical stresses at the air-wall interface of the human nasal cavity during breathing. J. Appl. Physiol. 100, 1003–1010. https://doi.org/10.1152/japplphysiol.01049.2005. Gambaruto, A.M., Doorly, D.J., Yamaguchi, T., 2010. Wall shear stress and near-wall convective transport: comparisons with vascular remodelling in a peripheral graft anastomosis. J. Comput. Phys. 229, 5339–5356. https://doi.org/10.1016/j.jcp.2010. 03.029. Girardin, M., Bilgen, E., Arbour, P., 1983. Experimental study of velocity fields in a human nasal fossa by laser anemometry. Ann. Otol. Rhinol. Laryngol. 92, 231–236. https://doi.org/10.1177/000348948309200304. Goodale, J.L., 1896. An experimental study of the respiratory functions of the nose. Bost. Med. Surg. J. 135, 457–460. https://doi.org/10.1056/NEJM189611051351901. Grgic, B., Finlay, W., Burnell, P.K., Heenan, a., 2004. In vitro intersubject and intrasubject deposition measurements in realistic mouth–throat geometries. J. Aerosol Sci. 35, 1025–1040. https://doi.org/10.1016/j.jaerosci.2004.03.003. Hahn, I., Scherer, P.W., Mozell, M.M., 1993. Velocity profiles measured for airflow through a large-scale model of the human nasal cavity. J. Appl. Physiol. 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. https://doi.org/10.1007/s00348-003-0636-7. Holbrook, L.T., Worth Longest, P., 2013. Validating CFD predictions of highly localized aerosol deposition in airway models: in vitro data and effects of surface properties. J. Aerosol Sci. https://doi.org/10.1016/j.jaerosci.2013.01.008. Hopkins, L.M., Kelly, J.T., Wexler, a.S., Prasad, a.K., 2000. Particle image velocimetry measurements in complex geometries. Exp. Fluids 29, 91–95. https://doi.org/10. 1007/s003480050430. Hornung, D.E., Leopold, Da, Youngentob, S.L., Sheehe, P.R., Gagne, G.M., Thomas, F.D., Mozell, M.M., 1987. Airflow patterns in a human nasal model. Arch. Otolaryngol. Head Neck Surg. 113, 169–172. Inthavong, K., Shang, Y., Tu, J., 2014. Surface mapping for visualization of wall stresses during inhalation in a human nasal cavity. Respir. Physiol. Neurobiol. 190, 54–61. https://doi.org/10.1016/j.resp.2013.09.004. Inthavong, K., Wen, J., Tu, J., Tian, Z., 2009. From CT scans to CFD modelling – fluid and heat transfer in a realistic human nasal cavity. Eng. Appl. Comput. Fluid Mechan. 3, 321–335. Ito, K., Mitsumune, K., Kuga, K., Phuong, N.L., Tani, K., Inthavong, K., 2017. Prediction of convective heat transfer coefficients for the upper respiratory tracts of rat, dog, monkey, and humans. Indoor Built Environ. 26, 828–840. https://doi.org/10.1177/ 1420326X16662111. 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. https://doi.org/10.1016/j.jaerosci.2008. 06.002. Kabilan, S., Suffield, S.R., Recknagle, K.P., Jacob, R.E., Einstein, D.R., Kuprat, A.P., Carson, J.P., Colby, S.M., Saunders, J.H., Hines, S.A., Teeguarden, J.G., Straub, T.M., Moe, M., Taft, S.C., Corley, R.A., 2016. Computational fluid dynamics modeling of Bacillus anthracis spore deposition in rabbit and human respiratory airways. J. Aerosol Sci. 99, 64–77. https://doi.org/10.1016/j.jaerosci.2016.01.011.

13