key: cord-0072092-baoy2fdf
authors: Wang, Li; Ge, Haiwen; Chen, Liang; Hajipour, Alireza; Feng, Yaning; Cui, Xinguang
title: LES study on the impact of airway deformation on the airflow structures in the idealized mouth–throat model
date: 2021-12-23
journal: J Braz
DOI: 10.1007/s40430-021-03324-7
sha: 8322cfdc9b2c05a8613162326b171e34e63135eb
doc_id: 72092
cord_uid: baoy2fdf

To investigate the impacts of upper airway deformation on the airflow structures, the airflow fields in the trachea are simulated using three geometrical models considering three different levels of airway deformations. Structured grids are used to create the high-quality grids. Large eddy simulation with the Smagorinsky sub-grid model is adopted to solve the three-dimensional in-compressible Navier–Stokes equations using the solver pisoFoam in the open-source CFD software OpenFOAM. The numerical results demonstrate that the airway deformation influences the main airflow structures depending on the deformation level. Particularly, it slightly impacts on the laryngeal jet such as the profile and the strength of laryngeal jet. The strength of the laryngeal jet increases slightly for the heavy deformation. In contrast, it impacts on the recirculation zone, secondary vortices, and turbulent kinetic energy more obviously. The increasing airway deformation will produce stronger secondary flow, smaller recirculation zone, and weaker turbulent kinetic energy. The turbulence intensity distribution varies as well. The obviously impacted flow region is mainly within the region of one to six tracheal diameters downstream the glottis.

The properties of the airflow fields in the human respiratory system are very complex and also very important to human health. For instance, coronaviruses that are delivered directly by the human respiratory airflow can cause serious threat to the public health system [1] . Particularly, COVID-19 has infected more than 0.17 billion people and has led to more than 3.8 million deaths globally until June 18, 2021 [2] . This has significantly increased the researchers' interests on the human respiratory multiphase flows [3] . Moreover, inhaled aerosol drug delivery via the oral/nasal airway also significantly benefits the human health to treat different kinds of diseases [4] . Therefore, well understanding airflow dynamics in the respiratory airways is very important since it is the foundation to understand the particle dynamics, which directly links to the drug delivery efficiency and the health damage of air pollution [4] .

There are mainly two kinds of geometrical models: the idealized geometrical models [5] and image-based patientspecific geometrical models [6, 7] . Although the realistic models based on the medical scans are adopted more and more recently due to the improvement in high-resolution computed tomography (CT) and geometry reconstruction technique [8] , the idealized models are still being used in this community. On one side, the complexity of the respiratory airways makes it impossible to resolve all of the bronchial airways, particularly for the small bronchioles [9, 10] ; on the other side, the idealized model is adopted when it focuses on studying the physical properties of the airflow dynamic since such kind of models has simpler airway wall boundary and retains the primary properties of the respiratory airways.

The airflow structures in such kind of cast-based mouth-throat airway model have been studied extensively. It is well known that the main airflow structures in the throat include the counter-rotating secondary vortices, recirculation zone, laryngeal jet, and vortex [4, 5, 11] . The laryngeal jet is very unsteady, and the tail of laryngeal jet breaks up [12] . The recirculation zone is highly unsteady as well, and it interacts with the laryngeal jet. Different length-scale secondary vortices distribute widely in the throat region. Relatively, the turbulent kinetic energy and turbulence intensity are not discussed so much. And few research has compared the properties of these flow structures when it considers different levels of airway deformation after glottis in this idealized cast-based model.

Usually, the wall of the respiratory airway is assumed to be solid without considering the airway deformation. However, the airway deformation or modification can occur during the process of breathing or in the special scenarios. Some investigations have reported the impacts of airway profiles on the airflow structures. For instance, Choi and Wroblewski [13] report that the airflow patterns are different with triangular glottis and circular glottis. Brouns et al. [14] also find that different shapes and cross section areas will lead to different airflow fields using RANS simulations. Xi et al. [15] report that the laryngeal jet instability and vortex generation vary when the glottal aperture deforms by varying the main flow speed. Paz et al. [16] study the cough clearance process by modeling the throat deformation based on a cast-based mouth throat model. Zhao et al. [17] compare inhaled drug deposition fraction of steady glottis and dynamic glottis. The glottis without deformation significantly deviates the total deposition fraction predictions. Zobaer et al. [18] simulate the airflow during the inhalation using trachea models with different degrees of stenosis. The results show that stenosed trachea will increase the required power for the patient to breathe. Jayaraju et al. [19] analyze the airflow in four different narrowing airways and find that the pressure drop increases when the degree of constriction increases. However, few research has been conducted to study the impact of the upper airway deformation after glottis on the airflow structures in detail. Xi et al. [20] generated infinitely large numbers of parameterized lung models with varying levels of airway distensibility (compliance) and constriction (resistance). The results show that local constriction can strikingly alter the airflow and particle deposition mapping throughout the lungs.

Numerical modeling is a powerful tool to analyze the airflow field, which provides more details of the airflow structures and has become a very popular way in this research area [4, 11, 12, 21, 22] . Nowadays, it is impossible to simulate the airflow field in the entire respiratory system due to its complexity of the respiratory system and the limitations of computational resources [9, 23] . The large difference in size through the whole respiratory system leads to a wide range of length scales. Therefore, most of the research studies are carried out in the different segmented airways such as the upper airways [24] [25] [26] [27] and tracheobronchial airways [11, 28, 29] . Due to its acceptable computational time and relatively low requirements of computational resources, the Reynolds-averaged Navier-Stokes equations (RANS) coupled with the low-Reynolds number k − turbulence model are the most popular way to predict the laminar-transitional-turbulent flows in the respiratory acts [30] presently. However, usually the RANS model is limited to predict the steady airflow structures. In contrast, the large eddy simulation (LES) [31] is capable to predict the unsteady flow structures and has been reported to better predict the transitional flow than the RANS model. For instant, Cui and Gutheil [12] have validated that LES can better capture the properties of airflow field via comparing the modeling results from RANS and LES with the experimental data. Therefore, although this turbulence model requires more computational resources and computational time, the LES method becomes more and more popular nowadays. Following the previous work [12, 25] , the LES is employed in the present work to predict the turbulent flows in the mouth-throat models. Although a lot of researchers adopt commercial software such as Ansys Fluent [30] to solve the governing equations of multiphase flows, OpenFOAM has its advantages such as feasibility and free access and thus has been adopted by several publications [25] In the present work, to study how the airway deformation can impact the properties of airflow field in the upper airway, the transient airflow field is simulated. Three geometrical models including one baseline model, one light throat deformation model, and one heavy throat deformation model are adopted to describe different levels of airway deformations. An idealized cast-based mouth-throat model is adopted, and the structured grid is employed to create the high-quality grids. The deformation of the airway after the glottis is given according to the work [16] . The complex laminar-transitional-turbulent airflow fields are modeled using a LES method. Smagorinsky sub-grid-scale model is adopted to solve the three-dimensional in-compressible Navier-Stokes equations using the solver pisoFoam in the open-source CFD software OpenFOAM. The paper will focus on discussing how airway deformations after glottis in the upper airway models impact the properties of the airflow fields in the trachea.

In this work, it assumes the airflow in the mouth-throat airway is incompressible due to its slow speed compared to the speed of sound. It also neglects the heat transfer and does not consider the impact of the temperature. Solid wall is adopted by excluding the effect of mucus layer. In Sect. 2, it will introduce the methodology including geometrical models, governing equations and boundary conditions, and numerical solution procedure. In Sect. 3, it will present and discuss the airflow properties under different levels of airway deformations such as the properties of laryngeal jet, recirculation zone, secondary flow, turbulent kinetic energy, turbulence intensity, and vortex. In the last section, a conclusion is made on this work.

An idealized mouth-throat model is created based on the human cast [5, 22, 24, 32] . This idealized mouth-throat model comprises the oral cavity, pharynx, larynx, and trachea. It has been widely adopted in the investigations of the multiphase flows in the human upper airway [4, 5, 22] . This model assumes that the cross section of the trachea is circular. However, as it is reported by Malvè et al. [33] , the posterior wall of the trachea is close to a flat plane due to the existence of cartilage rings [16] . Thus, we adopt the same way as their work to modify the posterior wall of the original mouth-throat model to make it closer to the reality, as shown in Fig. 1 . This modified model is set as the baseline model. In the real situations, high velocity due to coughing or other clinical complications that alter the trachea shape have been reported [33] [34] [35] . The light and heavy airway deformations after glottis are derived from the aforementioned baseline model as reported in the literature [16] , in which the posterior wall is deformed toward inside of the airway (c.f., Fig. 1 ). The deformed wall is half radius away from the central line for the case of light deformation and 0.25 radius for the case of heavy deformation. With this way, the impact of airway deformation on the properties of airflow fields after the glottis can be evaluated. The aforementioned three geometrical models, shown in Fig. 1 , are discretized with hexahedron cells using a commercial preprocessing software, Ansys ICEM-CFD. The grid independence was conducted previously by comparing the velocity along the central line in the similar geometrical model [12] . The similar strategy to generate the grid is adopted, and the final grid nodes are 1,276,500.

The properties of airflow fields in the upper airway are governed by three-dimensional in-compressible Navier-Stokes equations. To solve the transitional flow, LES is adopted and the sub-grid-scale (SGS) stress is modeled by the Smagorinsky model [36] . The Smagorinsky SGS model is easy to be implemented and widely used. Other SGS models will be implemented in the future. The filtered continuity and momentum equations are presented as follows: Continuity equation [31] Momentum equation [31] where Ū i presents the velocity component in the i direction and the bar presents the filtered variables; is the gas density; is the kinematic gas viscosity; p is the static pressure; and r ij is the sub-grid-scale stress tensor. In the Smagorinsky model [36] , the sub-grid stress tensor is modeled by where S ij is the strain rate tensor [36] Moreover, and the filter width is calculated from local cell size by = ( x y z) 1 3 in implicit filtering. C s is set as 0.167. In these equations, i, j = 1, 2, 3 and the Einstein summation convention is used.

Boundary conditions are set in the view of representing a typical inspiration process: The mouth inlet is prescribed by a fluctuating velocity with a mean flow rate of 30 L/min, which is corresponding to breathing intensity at light activity [11, 37] . No-slip wall boundary condition is applied to the airway wall; Jourabian et al. have validated that wall

function for large eddy simulation will cause simulation errors [38, 39] . A resolved no-slip wall boundary is applied to the airway wall. Outlet is modeled as an open boundary with zero gauge pressure. The density and kinematic viscosity of the air are = 1.21kg∕m 3 and = 15 × 10 −6 m 2 ∕s , respectively. In this work, it does not consider the heat transfer or the effect of the mucus layer on the investigation since it focuses on the airflow structures.

In this study, the governing equations are solved using the open-source software platform OpenFOAM [12, 25, [40] [41] [42] . Time-dependent term is discretized using the first-order implicit Euler scheme. Time step is set as 6 × 10 −6 second. The convective term is discretized by the limited secondorder central scheme that limits the flux toward upwind side in the regions of rapidly changing gradient. With respect to the diffusive scheme, a second-order central difference scheme with non-orthogonality correction components is employed. The pressure-velocity coupling is achieved using a modified pressure implicit with splitting of operators (PISO) scheme [43] . The modified Rhie-Chow [44] interpolation scheme is employed to construct a pseudo-momentum at the element face that avoids the checkerboard problem of the collocated grid. All of the numerical cases are run at the High-Performance Computing Center (HPCC) in the Texas Tech University. The computers are Dell PowerEdge. The processors are the Intel Xeon E5. Each core has a memory of 5.33GB. For each case, 180 cores are used. The maximum Courant number is maintained to be below 0.99, and it takes about 350 hours to run the simulation of flow time 3 s.

In this section, airflow structures are discussed for the breathing rate at 30 L/min. The corresponding Reynolds number is 2,320 at the inlet of the mouth cavity. The dimensionless velocities are adopted to compare the properties of the main airflow structures with different airway deformations, which are normalized by the velocity magnitude at the mouth inlet. Figure 2 displays the contours of the velocity magnitude at the mid-planes of the idealized mouth-throat models. Figure 3 displays the magnitude of the axial velocity contour at different cross sections. The four cross sections A-D refer to positions at the glottis, one, three, and six diameters of the trachea downstream the glottis that are indicated in Fig. 2 . Figure 4 displays the velocity magnitude and streamlines of the secondary flow at the cross sections A-D. Figure 6 displays the three-dimensional profiles of laryngeal jets using ISO surface of the total velocity magnitude. Figure 7 displays the recirculation zones using the three-dimensional flow streamlines marked with velocity in the y direction. In these figures, U, U y , and U s indicate the averaged velocity magnitude, averaged axial velocity magnitude, and averaged velocity magnitude of the secondary flow, respectively. The letters 'U' and 'u' indicate averaged and instantaneous airflow field, respectively. From these figures, it can be seen that the present simulation can capture the main flow structures including secondary vortices, laryngeal jet, and recirculation zone in the upper airway as reported in [5, 12, 22] . Moreover, both these flow structure profiles and quantities are very close to the results in [5, 22] from another research group, which can sufficiently validate the present modeling. In addition, the present methodology is further validated in [12] by comparing the numerical result with the experimental result in a constricted tube. The authors notice that it requires more rigorous and relevant validation with the experimental data in the future (Fig. 5 ).

As pointed out by the previous studies [12, 19] , the narrower glottis will induce a laryngeal jet with higher velocity due to mass conservation, which is shown in Fig. 2 . From the plots in the top row of this figure, it can be seen that the profile of the laryngeal jet is very similar in the averaged airflow field for all the cases. (Hereafter, it refers to baseline case, light deformation, and heavy deformation, respectively, if it only mentions all the cases.) The highest velocity levels are 7 (dimensionless velocity, hereafter) for all the cases. The profiles of the contour of axial velocity magnitude at cross sections A and B in Fig. 3 (top) further confirm that the strengths and profiles of laryngeal jets are similar for all the cases. All of the highest axial velocity magnitudes are at the levels of 7 at the cross section A and 6.5 at the cross section B. The deformation has negligible impacts on the flow field at the upstream of cross section B. When the airflow goes further downstream to the cross section C that is downstream of the deformation location, the impacts of the airway deformation become evident: The dimensionless axial velocity is 6.5 for the case of heavy deformation, and 6 for the baseline and light deformation cases. This is because the heavy airway deformation reduces the cross area of throat more and leads to larger velocity due to the mass conservation law. At cross section D, the axial velocity becomes much more homogeneous after the turbulent flow In the instantaneous airflow fields (c.f. Fig. 2 (bottom) ), the profiles of laryngeal jets are very similar before cross section B. However, they are evidently different in the tails of the laryngeal jets. From Fig. 3 (bottom) , it can further be observed that the profiles of laryngeal jets are similar for all the cases at cross sections A and B, while much more different at cross sections C and D. This further confirms that the deformation has negligible impacts on its upstream airflow field. In addition, as same as the averaged airflow field, the instantaneous velocity is higher with increasing throat deformation, whose highest levels are 3.5, 3.8, and 4.0 for the baseline, light deformation, and heavy deformation cases, respectively.

To further investigate the laryngeal jet, ISO surfaces of the velocity magnitude with values of 4 and 4.5 are plotted for the averaged and instantaneous airflow fields, respectively (c.f., Fig. 6 ). All of the laryngeal jets are concave, and the profiles are very similar in the averaged airflow fields, while the laryngeal jet is more concave for the case of heavy deformation than the baseline case (c. f. Fig. 6 (top) ). The length of these laryngeal jets is 4.12, 4.27, and 4.85 tracheal diameters, respectively, which are calculated from the glottis. Therefore, the length of the laryngeal jet slightly increases when the airway deformation becomes heavier. From Fig. 6 (bottom) , it also shows that the laryngeal jets are very unsteady and the tails of laryngeal jets are breaking up. Moreover, from this figure, we can further confirm that the laryngeal jet is much steadier for the heavy deformation case in the upstream than the other cases.

In this section, another important airflow structure-the recirculation zone-in the trachea is also discussed. As shown in Fig. 2 , the lower velocity zones in the posterior side of trachea display the recirculation zones. The recirculation zones are generated due to the expansion of the cross area in the larynx and throat. This expansion will generate adverse pressure gradient, and then the recirculation zone forms in the throat. The airway deformation will reduce the cross area and then reduce the recirculation zone as seen in cross sections B and C in Figs. 2 and 3. It can be clearly observed from Fig. 2 that the mixing between the recirculation zone and the main flow region becomes smaller when the airway deformation increases. From Fig. 7 , it can be seen that the velocity magnitude of the recirculation zone is lower when the airway deformation is heavier. This figure also demonstrates that the recirculation zone locates at the posterior side of the throat for the baseline case, while it shifts to the corners when the airway deforms. These indicate that the airway deformation will modify the location and profile of the recirculation zone. In addition, the streamlines Fig. 7 ). It indicates that the recirculation zone is much more unsteady for the baseline case than the case of heavy deformation. When the deformation is heavier, the characteristic length scale becomes smaller. Thus, the corresponding Reynolds number becomes smaller and the flow becomes less turbulent.

In summary, the airway deformation has a slight impact on the strength and profile of the laryngeal jet. The impact is much more obvious at the downstream of the laryngeal jet and for the case of heavy airway deformation. The size of recirculation zone is reduced, and the profile is modified when the airway deforms.

The secondary flow is important since it affects the particle transport and deposition. From Fig. 4 , it can be seen that there is a pair of counter-rotating vortices in the averaged airflow field locating at the anterior side of glottis at the cross section A for all the cases. The lower velocity magnitude of the secondary flow locates at the region of secondary vortices core. The higher velocity magnitude of secondary flow locates at the middle region of the cross section due to larger shear stress in this area. There are no evident differences at this section among all the three cases in the terms of both profiles of the secondary vortices and the magnitudes of the secondary flow. When the airflow goes further downstream, the secondary flow becomes a little more complex with another pair of secondary vortices appearance at the cross section B for the baseline case and the heavy deformation case, and another two pairs of the secondary vortices for the light deformation case. Meanwhile, the velocity magnitudes of the secondary flow are very different near the boundary of the recirculation zone and the laryngeal jet at the cross section B. While the highest levels of secondary flow velocity magnitude are still 2 at cross section B, the region of high velocity magnitude becomes smaller in comparison with the one at cross section A. At the cross section C, a pair of large-size secondary vortices exists at Fig. 6 Iso-surface of the velocity magnitude presenting the laryngeal jets in the averaged (top) and instantaneous (bottom) airflow fields for the cases of base, light deformation, and heavy deformation this section for the baseline case. In contrast, in addition to the large-size counter-rotating vortices, a pair of small-size vortices is observed at this section for the cases of light and heavy airway deformations. Coupling with differences of the secondary vortices, secondary flow velocity magnitudes are different as well. The obvious differences are demonstrated clearly in the posterior side of the trachea as well as in the vortex core region shown in Fig. 4 (top) . The strength of the secondary flow decreases at cross section C in comparison with the flow at cross sections A and B. The highest velocity magnitude of secondary flow at cross section C is larger for the heavy deformation than the ones of the light deformation and baseline cases due to the reduction in the cross section area. At cross section D, the secondary flow becomes much weaker for all the cases. The magnitude of the secondary flow is 0.2, 0.3, and 0.35, respectively. At this cross section, the differences of the secondary vortices among the three cases are more evident. Although velocity distributions of secondary flows are similar at this cross section for all the cases, the values and profiles are evidently different, particularly at the posterior and lateral sides of the throat. The size of the smaller vortices increases for both the light and heavy deformation cases at the cross section D in comparison with cross section C.

In the instantaneous airflow field (c.f. Fig. 4 (bottom) , the highest levels of the secondary airflow velocity magnitude are 2 at the cross section A, which are the same as the averaged airflow field. The profiles of the secondary vortices are different between the averaged and the instantaneous airflow fields, while the large-size counter-rotating vortices are still be observed. The similar velocity contours are observed between the averaged and instantaneous flow field at this cross section for all the cases. This indicates that the flow is less unsteady at glottis (cross section A) than the flow in the throat seen at cross sections B-D. It implies that the airflow properties at cross section A are not significantly influenced by the airway deformation. At cross sections B and C, the profiles of the secondary vortices are very different from the ones of averaged flow field for all the cases. Also, the instantaneous velocity magnitude contours of secondary flows are different from the averaged ones for all the cases. At the cross section D, the secondary vortices randomly distribute at this section. It can be seen that the highest levels of the instantaneous secondary flow velocity magnitudes are 2.4, 2.0, and 3.0 at the cross section C and 0.9, 1.0, and 1.3 at the Fig. 7 Recirculation zones are demonstrated via the contoured flood of the streamlines with the axial velocity in the averaged (top) and instantaneous (bottom) airflow fields for the cases of base, light deformation, and heavy deformation cross section D for all three cases, respectively. In comparison with the ones of the averaged airflow fields, the highest levels of the secondary velocity magnitude are very high for the instantaneous airflow field. This indicates that the fluctuations of secondary flow are very large.

To further analyze the properties of the secondary flow, the velocity magnitude of secondary flow and their ratio to the total velocity magnitude is illustrated in Fig. 5 cross section) , respectively, where d is the tracheal diameter. For the secondary flow along the line 1, the velocity magnitude of the secondary flow arrives peak near the location of 0.52 (dimensionless distance with tracheal diameter to the glottis, hereafter) and then decreases for all the three cases. After the location of 1.3, the secondary velocity and its ratio are larger for the heavy deformation case than the other two cases. There is a lowest region at the distance of 1.9 for line 2 and 1.8 for line 3 for both the velocity magnitude of the secondary flow and its ratio to the total velocity magnitude, as shown in Fig. 5 . At the line 3, the profile of the secondary flow velocity magnitude and its ratio to the total velocity magnitude are obviously different for the heavy deformation case from the other two cases, particularly at the locations between 1 and 4. After 4, both of the velocity magnitude of the secondary flow and its ratio are larger for the heavy deformation case than the ones of the other two cases. It indicates that the airway deformation will alter the distribution of secondary flow and will generate more mass flux in the secondary flow that leads to more swirling flow. In summary, the airway deformation has a significant impact on the airflow field, which has been demonstrated by the secondary vortices distribution, velocity magnitude of secondary flow, and its ratio to the total velocity, in particular for the heavy airway deformation.

The dimensionless turbulent kinetic energy k is turbulent kinetic energy ( 1 2 (u �2 + v �2 + w �2 ) [45] ) normalized by U 2 I , where u ′2 ,v ′2 ,w ′2 are root-mean-square components of local velocity and U I is the inlet velocity, dimensionless turbulent kinetic energy. Their contours are displayed at the mid-plane and cross sections A-D as shown in Fig. 8 . From Fig. 8 (top), it can be seen that the highest levels of turbulent kinetic energy ( k ) at the mid-plane are 1.8, 1.6, and 1.6 for all the cases, respectively. It means that the increasing airway deformation will decrease the turbulent kinetic energy, which is consistent with the properties of the recirculation zones. In the region between A and B cross sections, the turbulent kinetic energy is lower for the case of the heavy airway deformation, which is consistent with previous observations that the airflow field is much steadier in this region for the case of heavy deformation. The turbulent kinetic energy is mainly produced in the regions of tracheal airway between cross sections B and D. As shown in Fig. 8 (bottom), much smaller k is seen at the cross section A than the other three cross sections. The profile of k at cross section A is almost the same for all three cases, whose highest levels are 0.25. At the cross section B, evidently different profiles of k are observed. The highest levels of k at this cross section are 1.8, 1.6, and 1.0, respectively. The larger k appears at the boundary of the laryngeal jet, where the shear stress is larger and the production of turbulent kinetic energy is higher. The maximum k locates in the posterior side of the throat for the baseline and light deformation cases. The maximum of k locates in the lateral sides of throat for the heavy deformation case. This indicates that the airway deformation has modified the way of the production of k. At cross section C, the highest levels of k are 1.5, 1.3, and 1.4 for the three cases, respectively. The larger k locates at the boundary of the recirculation zone and its interface with the laryngeal jet at the cross section C. The distributions of k at this section are obviously different from the other cross sections as shown in Fig. 8 (bottom) . At the cross section D, the highest levels of k decrease to 0.36, 0.3, and 0.36 for the three cases, respectively. k distributes evenly at this cross section (c. f. Fig. 8 (bottom) ). To further investigate the variation of k, k along three axial lines is plotted in Fig. 10 , in which the legends and the locations of the axial lines are the same as Fig. 5 . It can be seen from these plots that the k increases at first and then decreases for all the lines and all the cases. The turbulent kinetic energy is higher for the baseline and light deformation cases in the most regions, in particular, in the regions within one and four tracheal diameters downstream the glottis (c.f., line 2 and line 3). However, the k is higher near the peak for the case of heavy deformation than light deformation at line 1. From Fig. 10 , it can also observe that the k is higher in the lateral position (line 3).

Furthermore, the turbulent intensities ( TI = √ 2∕3kU 2 I U ) [45] are stronger in the recirculation zones and their interface with the laryngeal jet seen at the mid-plane (c.f. Fig. 9 (top)) and at the cross sections B and C (c. f. Fig. 9 (bottom)). The turbulent intensities at cross sections A and D are relatively smaller than the other two cross sections. The distributions of turbulent intensities at all cross sections are similar among different cases. It can be seen that the TI distributions are not always consistent with the turbulent kinetic energy distributions, i.e., at the cross sections B and C (c. f. Fig. 9 (bottom) ). According to Fig. 10 (bottom) , the peak locations of the TI are different from the ones of k for all the cases. This is because TI is not only related to k distribution, but also with the averaged velocity. Generally speaking, the TI of the baseline case is larger than the ones of airway deformation cases. However, the TI is larger for the heavy deformation case in the peak area (seeing line 1 and line 2 in Fig. 10) .

In this study, the third-generation vortex identification method called the -criterion method [46] is adopted to investigate the vortex. Compared with the previous studies using the Q-criterion or the -criterion, it defines the vortex concentration [46] as where is a small number to avoid division by zero, and [46] and [46] are the symmetric part and anti-symmetric part of the velocity gradient tensor ∇ , respectively: This method has advantages to capture more realistic vortex. Figure 11 demonstrates the 3D vortex structures in the instantaneous airflow fields for the three cases. The 3D vortex structures are identified by -criterion with =0.6 for all the cases. There are more large-size vortices in the cases of light and heavy deformations than in the baseline case, particularly in the region between the glottis and B cross section. These are consistent with the observations about

velocity contour, secondary streamlines, and the turbulence intensity. The higher k means the airflow is more unsteady. The higher turbulence intensity in the baseline case than the one in the heavy deformation case indicates that the velocity fluctuation is larger. The higher turbulence intensity implies that more kinetic energy is converted into the turbulent kinetic energy that eventually will be dissipated into thermal energy. In summary, the airway deformation after glottis in the upper airway will impact the properties of the turbulent kinetic energy, the turbulence intensity, and the turbulent vortex in the tracheal Three-dimensional turbulent vortices displaying with = 0.6 for the cases of base, light deformation, and heavy deformations airway, which is mainly because it has modified the recirculation zone and the shape of the laryngeal jet.

In this study, the impacts of airway deformation on the properties of airflow field in the human respiratory airway are investigated using a LES method with the airflow rate of 30 L/min (Re = 2,320) at the inlet of the mouth cavity. The properties of the main airflow structures such as the laryngeal jet, the recirculation zone, and the secondary flow are analyzed by considering different levels of airway deformations. The turbulent kinetic energy, turbulence intensity, and vortex are discussed as well. The following conclusions can be drawn from the present study:

(1) The airway deformation after glottis in the upper airway has minor impact on the laryngeal jet if the deformation is light. However, it will slightly impact the profile and length of the laryngeal jet if the deformation increases.

The increasing deformation will reduce the size of the recirculation zone, and the interaction of the laryngeal jet and recirculation zone will become weaker. (2) The properties of the secondary flow are closely related to the airway profile. The airway deformation will modify the profiles of the secondary vortices as well as the strength of secondary flow. With the airway deformation increasing, the secondary flow becomes stronger. (3) The turbulent kinetic energy becomes weaker when the airway deformation increases. Due to the variation of the recirculation zone and the interface between the recirculation zone and laryngeal jet, the properties of k as well as the turbulence intensity are modified. Meanwhile, more three-dimensional large-size vortices are observed in the heavy deformation case.

To better capture the airflow field, more realistic geometrical model based on CT scans, heat transfer as well as effects of mucus layer will be adopted. Additionally, more advanced SGS model or even DNS method will improve the accuracy of the numerical simulations. These will be considered in our future research work.

Update on human coronaviruses: one health, one world

Coronavirus disease (covid-19) pandemic

Numerical modeling of the distribution of virus carrying saliva droplets during sneeze and cough

Targeted drug-aerosol delivery in the human respiratory system

Computationally efficient analysis of particle transport and deposition in a human wholelung-airway model. part i: Theory and model validation

Acoustic simulation of a patient's obstructed airway

Numerical simulation of the dispersion and deposition of a spray carried by a pulsating airflow in a patient-specific human nasal cavity

A review of respiratory anatomical development, air flow characterization and particle deposition

Computational analysis of aerosol-dynamics in a human whole-lung airway model

Use of computational fluid dynamics deposition modeling in respiratory drug delivery

Airflow and particle transport in the human respiratory system

Large eddy simulation of the unsteady flow-field in an idealized human mouth-throat configuration

Characteristics of glottis-induced turbulence in oscillatory flow: an empirical investigation

Influence of glottic aperture on the tracheal flow

Effects of glottis motion on airflow and energy expenditure in a human upper airway model

Cfd transient simulation of the cough clearance process using an eulerian wall film model

Glottis motion effects on the particle transport and deposition in a subject-specific mouth-totrachea model: a cfpd study

Modeling the effect of tumor compression on airflow dynamics in trachea using contact simulation and cfd analysis

Effects of tracheal stenosis on flow dynamics in upper human airways

The application of statistical shape modeling for lung morphology in aerosol inhalation dosimetry

Planning human upper airway surgery using computational fluid dynamics

Fluids moving boundary simulation of airflow and micro-particle deposition in the human extra-thoracic airway under steady inspiration. part i : airflow

Particle deposition in a realistic geometry of the human conducting airways: Effects of inlet velocity profile, inhalation flowrate and electrostatic charge

Laminar-to-turbulent fluidparticle flows in a human airway model

Three-dimensional unsteady large eddy simulation of the vortex structures and the mono-disperse particle dispersion in the idealized human upper respiratory system

Large eddy simulation of polydisperse particle deposition in an idealized mouth-throat

Numerical study of the airflow structures in an idealized mouth-throat under light and heavy breathing intensities using large eddy simulation

Characteristics of the turbulent laryngeal jet and its effect on airflow in the human intra-thoracic airways

Direct numerical simulation of particle laden flow in a human airway bifurcation model

Effects of generation time on spray aerosol transport and deposition in models of the mouth-throat geometry

Large eddy and detached eddy simulations of fluid flow and particle deposition in a human mouth-throat

Particle deposition in a cast of human oral airways

Fsi analysis of the coughing mechanism in a human trachea

Modeling of the fluid structure interaction of a human trachea under different ventilation conditions

Elastic properties of the central airways in obstructive lung diseases measured using anatomical optical coherence tomography

General circulation experiments with primitive equations-i, the basic experiment

Laminar-to-turbulent fluid-nanoparticle dynamics simulations: model comparisons and nanoparticle-deposition applications

Large eddy simulation (les) of suspended sediment transport (sst) at a laboratory scale

Wall-layer model for large eddy simulation (les) of suspended sediment transport (sst) in a labscale turbulent open channel flow

Numerical and experimental study on the deposition of nanoparticles in an extrathoracic oral airway model

Laminar/turbulent airflow and microsphere deposition in a patient-specific airway geometry using an open-source solver

reactingfoam-sci: An open source cfd platform for reacting flow simulation

Solution of the implicitly discretised fluid flow equations by operator-splitting

Numerical study of the turbulent flow past an airfoil with trailing edge separation

Study of the flow unsteadiness in the human airway using large eddy simulation

New normalized rortex/vortex identification method

Acknowledgements The authors acknowledge the High-Performance Computing Center (HPCC) at Texas Tech University at Lubbock for providing HPC resources that have contributed to the research results reported within this paper.

The authors state that there is no personal/commercial relationships involved in this work, which would lead to conflict of interest.