Turbulent Kinetic Energy of Flow During Inhale and Exhale to Characterize the Severity of Obstructive Sleep Apnea Patient

This paper aims to investigate and present the numerical investigation of airflow characteristics using Turbulent Kinetic Energy (TKE) to characterize the upper airway with obstructive sleep apnea (OSA) under inhale and exhale breathing conditions. The importance of TKE under both breathing conditions can show an accurate method in expressing the severity of flow in sleep disorder. Computational fluid dynamics are used to simulate the upper airway's airflow via steady-state Reynolds-averaged Navier-Stokes (RANS) with k–ω shear stress transport (SST) turbulence model. The three-dimensional (3D) airway model is created based on the CT scan images of an actual patient, meshed with 1.29 million elements using Materialise Interactive Medical Image Control System (MIMICS) and ANSYS software, respectively. Both high TKE and Rev were noticed around the region after the necking (smaller cross-sectional area) during the inhale and exhale breathing. The turbulent kinetic energy could be used as a valuable measure to identify the severity of OSA. This study is expected to provide a better understanding and clear visualization of the airflow characteristics during the inhale and exhale breathing in the upper airway of patients for the medical practitioners in the OSA research field.


Introduction
Obstruction sleep apnea (OSA) affects the quality of sleep for the patients. The OSA results in a temporary blockage of a patient's airway by several obstructions, such as collapse or laxity of airway tissues, which prevent normal airflow during breathing. The obstruction of the airway leads to hypopnea or apnea. Apnea is the complete blockage of air, while hypopnea is the partial blockage of air. Several situations differ the severity of blockage: relaxation of throat muscle, excessively bulky airway tissues, thick fatty tissue, and tongue falling back to the throat area. Apnea-Hypopnea Index (AHI) is practically used to determine the severity of OSA. Clinical data are collected to determine the severity of OSA via the average number of hypopnea or apnea events by every hour of sleep.
Previously, Faizal et al. [1] showed the importance of TKE in visualizing the accurate method in expressing flow concentration during light and heavy breathing for OSA patients.
The high intensity of TKE shown in flow concentration during light and heavy might induce high vibration to the soft tissue around the upper airways. Due to this, many previous researchers discovered that the vibration of soft tissue and snoring phenomena during a night of sleep might cause damage to the soft tissues, vessels, and nerves, thus resulting in sleep apnea for the patient [2,3]. However, the metrics of vibration and snoring are rarely correlated with the airflow effects. According to fluid mechanics principles, a turbulent airflow in a passage could cause vibration, leading to structural damage. This phenomenon has already been established in other engineering disciplines, such as in nuclear plants [4], aircraft design [7,8], and structural design [9]. Thus, it could be very well argued that a turbulent flow would cause vibration that could damage soft tissues, vessels, and nerves.
Computational fluid dynamics (CFD) has emerged as one of the powerful tools for attaining a deep insight into many medical disorders, including that of OSA [10]. Recently, many researchers focused on fluid mechanics in understanding the flow mechanism inside the upper airway [11]. It has been applied to study the upper airway anatomy and airflow and discover the virtual reality surgery in sleep-disordered breathing (SDB) [14,15]. Recently, Faizal et al. [16] reviewed the importance of computational simulation for understanding the human upper airway (HUA) in sleep apnea assessment. Powell et al. [17] constructed a complete CFD model that includes the details of the pharyngeal to investigate the human upper airway that has been used for surgery. Moyin et al. [18] utilized CFD to establish the pharyngeal airway model before and after surgery. The generalized airflow and pressure profile were used to determine the possible relationship with the treatment outcome. The research mentioned above considered seven OSA patients to investigate the upper airway. It was found that two OSA patients still experienced OSA even after treatment.
The flow pattern in the human upper airway [19] has been studied using experimental means or via CFD, where the results are focused on the velocity, pressure, and wall shear stress. Similarly, Mylavarapu et al. [20] used the CFD model to simultaneously model the airflow velocity, pressure, and wall shear stress in the human upper airway. With the knowledge of turbulent flow in the human upper airway, it is possible to avoid OSA. Therefore, the current study is the continuity research from our previous finding [1], and it is motivated to investigate the turbulent airflow in the pharyngeal airway for a patient diagnosed with OSA, mainly focusing on the inhale and exhale flow condition. The CFD simulation was applied for the task undertaken. The novelty is explained by visualizing the turbulent kinetic energy (TKE) to relate the severity of OSA patients with eddies of the airflow in the pharyngeal airway during the inhale and exhale breathing. The CFD result indicates a more reliable but relatively simplified upper airway model to understand OSA.

Method to Studying the Human Upper Airway
The method used in this work includes CFD modeling and segmentation of computed tomography scan (CT scan) acquisition been described and explained in previous research [1]. Therefore, Figure 1 portrays the CFD modeling in which the model was considered to segmentation the pharyngeal in the human upper airway. The volume of the model was created based on real patient computed tomography (CT) images. The focus was on both flow simulations. Thus, the marked locations of P6 and P1 are defined as the inlet and outlet boundary conditions for inhaling, but for exhaling, P6 is defined as outlet and P1 as an inlet. TKE in the upper airway was obtained in the simulation by considering all essential parameters and boundary conditions adopted in previous research.  The summarized characteristics of the OSA subject that were selected for the conduction of the upper airway CT scan. The focus of the subject of interest was decided to be from the nostril until the pharyngeal, according to previous works [18,21]. A total of 431 number of frames with 0.3 mm of slice thickness covering the upper airway CT scan of the subject were performed using an i-CAT Cone Beam 3D Dental Imaging System (version 3.1.62 supplied by Imaging Science International, Hatfield, USA), as shown in Figure 2. The airway boundary is identified from the CT scan images via a threshold based on the intensity of the gray image. The images were captured while the patient was awake and lying down with their face upward, as explained in the works of Chan et al. [22]. Each CT scan has 534 pixels ×

Mesh generation
The model generated by MIMICS software was transferred to 3-matic (version 15.0; Leuven, Belgium) to create the volume mesh. The preprocessor, ICEM (ANSYS, USA), generated unstructured tetrahedral meshes for airway volume. The mesh quality was controlled to ensure the highest level of result accuracy. The surfaces of the pharyngeal airway have meshed with 2.0 mm of Maximum face size. The maximum and minimum size of the mesh is 2 mm and 0.002 mm, respectively. The inflation for pharyngeal airway mesh was considered around the wall boundary to increase the accuracy of the simulation. This consideration also ensures an excellent fine mesh is created next to the wall, thus providing an adequately small mesh close to the size of the wall mesh. The y + (dimensionless distance) approach was used to accurately estimate the initial boundary layer thickness [17] that corresponds to the near-wall cell size (0.272). The quality of the meshed airway model was examined in ANSYS, where the quality of the mesh is measured by the average skewness value and standard deviation. The skewness value is 0.2142, which indicates an excellent cell quality (0-0.25), as recommended by ANSYS. Then, the mesh sensitivity grid study was performed on the 3D airway model with different grid scales. The 5-grid size was performed to establish the level of acceptable accuracy, as shown by the previous research approach.

Numerical modeling of the upper airway
Inhalation involves the airflow passing through the pharyngeal airway. The airflow characteristics may vary due to the uniqueness of the pharyngeal airway. Thus, in the simulation analysis, the flow field in the pharyngeal airway for OSA patients is computed by using a steady-state Reynolds-averaged Navier-Stokes (RANS) formulation, with the k-ω shear stress transport (SST) turbulence model with low Reynolds number correction [23].
The popular k-ω SST turbulence model was considered in this work instead of the k-ϵ model due to the greater accuracy of the former in viscous near-wall region treatment while considering the effects of adverse pressure gradient [24].

Governing equation
The continuity equation (equation 1) and momentum equation (equation 2) was considered and solved for the motion of incompressible airflow in the simulation analysis. However, the temperature change of the airway was minimal in this study. Thus, the energy equation was neglected.
Where, and are the density and viscosity of air, respectively; ⃗ is the air velocity vector; is the air pressure; and ⃗ is the gravity vector.

Turbulent kinetic energy (TKE) model and turbulent Reynolds number
The TKE equation is developed from the mass and momentum conservation equations to describe how the mean flow feeds kinetic energy into turbulence, given its essential role in developing the turbulent model in CFD modeling. The Reynolds decomposition method was used to decompose the velocity signal, , into mean ̅ and fluctuation component, ′ .
Therefore, the amount of turbulent fluctuation is described by the root-mean-square of the difference of instantaneous velocity and mean velocity [25], as represented below: where, is the number of samples in the signal included in the ANSYS simulation. The equation of fluctuation component is also applicable for the velocity signal, and . With these fluctuation components ( ′ , ′ , and ′ ), it is then possible to define the turbulent kinetic energy (TKE) as: where, ′ , ′, and ′ are fluctuating velocity components; and is the density. Using this equation, the velocity fluctuation used in calculating the TKE around the upper airway may be useful to identify the severity of blockage breathing in the upper airway. In comparison, kinetic energy (KE) is defined as: where, , , and are the phase-averaged velocity.
The definition of turbulent Reynolds number, [26] is shown in Eq. (6): The balance of kinetic energy is given in Eq. (5), while the work friction force, , is given in Eq. (7): where, is the total shear stress component that includes both molecular and turbulent friction to the surface of the integration volume, .

Boundary conditions
The boundary conditions of the 3D model were defined via the surface of the model, as mentioned in Section 2.0 (Figure 1). The inlet and outlet boundary conditions in the simulation were referred from the works of Mihai et al. [27]. The steady-state turbulent airway flow is computed by ANSYS Fluent, which solves finite volume using fitted grids.
During breathing, the velocity profile was assumed to be uniform, and the axial component of the velocity was perpendicular to the flow inlet, the nasopharynx [17]. Light breathing was considered at the inlet for both breathing conditions, where a volume flow rate of 7.5 L/min was considered in the simulation [28]. A turbulent intensity of 10% was considered to mimic the real condition, where the outlet was defined by an average gauge pressure of 0 Pa [29,30].

Grid sensitivity and solver verification
The grid sensitivity has been considered with a different total of elements to determine the sensitivity of the grid before proceeding in solver verification. After that, this research carries out the solver verification discussed in previous research [1]

Limitation of study
The methods were carried out in accordance with the approved guidelines of SCIENTIFIC REPORTS. Some of the limitations of this study include the limited extent of the airway understudy, the fact that the velocity of breathing is considered uniform, patients were awake and in an upright posture during the scans, the absence of a sleep study, the lack of sleepiness scales or questionnaires, the unknown stage of respiration, and the lack of tongue control. In addition, the airway is not a rigid structure, and the effects of the soft tissue placidness could not be taken into consideration.

Results and Discussion
The CFD simulation result of airway analysis is significant and valuable to a medical practitioner because it provides a clear airflow visualization. The current airway analysis assessment focuses on various aspects, such as the cross-sectional area, airflow patterns, velocity, pressure, wall shear stress, turbulent kinetic energy, and turbulent Reynolds number. Figure 3 shows the variation of airway cross-sectional area in thirty-four-plane images, with the location of the cross-sectional plane in the z-direction (i.e., flow direction). This indicates that the patient has a minimum cross-sectional area of 7.68 mm 2 at 28 mm from the inlet.

Cross-sectional area
The minimum area is crucial for patients with sleep apnea. Contraction of the crosssectional area influences the smoothness of the airflow when passing through the upper airway. This situation is also defined as the obstruction of airflow during inhalation. The narrow airway area in a patient with OSA might have been caused by the collapse or laxity of airway tissue. As mentioned beforehand, the RANS k-ω SST formulation was applied to simulate the airflow for a patient with OSA during inspiration of normal breathing at 30 L/min, a similarly used value by Powell et al. [17] to investigate the flow pattern in a pharyngeal airway. Most researchers, such as Zhu et al. [26], Hamidreza et al. [27], and Jernej et al. [28], have discussed and employed the flow values (i.e., velocity, pressure, and wall shear stress) to describe the characteristics of the upper airway. In this work, the additional flow parameter, TKE, is considered to describe the characteristics of the upper airway. Other flow parameters were also considered to visualize the phenomena and flow patterns in the upper airway. The results show that the highest velocity began after the airflow at the oropharynx. As the airflow approached the minimum cross-sectional area of the airway (retropalatal pharynx), the velocity increased, and the airflow seemed to have developed a turbulent flow pattern that also generated the flow recirculation regions (negative axial velocity), as obviously shown in Figure 4(b). The dilation of the cross-sectional area after the narrow area had affected the airflow in the OSA patient's airway. Thus, the flow recirculation is concentrated around this region, typically formed downstream with obstructions (narrow area). The existence of negative axial velocity may lead to the vibration of airway tissue due to continuous breathing.
Hence, it is highly likely that this region will generate vibration, causing the snoring. Snoring is known to reduce the strength of the tissue in patients with sleep apnea [3].   The TKE value will be high if the velocity fluctuation is in a randomized direction, with high magnitude impacting the blockage in the breathing area, as visualized in Figure 4; there is low pressure at the narrow cross-sectional breathing area, especially during the exhale. The high wall shear stress is mainly observed in the obstructed regions, contributing to high turbulent energy near the wall [25]. In a recent study, Xi et al. [31] found that shear stress generates fluctuating velocity and produces TKE in this region. The fluctuating velocity and turbulent flow near the wall will induce vibration and simultaneously generate snoring. As such, the severity of snoring directly corresponds to the TKE. The occurrence of snoring in the long term would cause reduced tissue strength for the pharyngeal, which will lead to the laxity or collapse of tissue due to gravitational effects when the patient sleeps in a lying down position. Consequently, this makes it a contributing factor to sleep apnea [2].

Turbulent kinetic energy (TKE) distribution along the pharyngeal airway
The turbulent kinetic energy is associated with the kinetic energy of vortices (velocity fluctuation) generated in the turbulent flow, which quantifies the turbulence level in the airway. In the clinical study, the index to measure the severity of OSA is the apnea-hypopnea index called AHI, a diagnosis of the periodic breathing of a patient during sleep. Hypopnea is when the patient undergoes shallow breathing due to constricted airway. The constricted airway will induce the intensity of the flow, causing high turbulent flow that creates velocity fluctuation, as shown in Figure 9. The fluctuation of velocity (Eq. (3)) and TKE (Eq. (4)) are related to hypopnea, which causes shallow breathing during sleep in measure of the AHI and TKE in CFD simulation.
The CFD analysis computes TKE and turbulent Reynolds number Rey. Figures 7 and 9 show the distribution of TKE obtained from the velocity fluctuation, as shown in Eq. (4); whereas the turbulent Reynolds number was computed from Eq. (6), following the works by Vladimir et al. [26] and Mohit et al. [35]. The airflow motion in the airway is predominantly by the energy in turbulent flow; therefore, the application of TKE is deemed suitable to describe the vortex and quantify the turbulence level. The results show high turbulent kinetic energy concentrated after the narrow region. In addition, the formation of a jet stream was observed in the pharynx of the airway. The jet stream formed had hit the wall in the inhalation area, which may have caused velocity remodeling near the airway wall. The advent of high-resolution three-dimensional TKE able to be constructed from a synthesis of standard two-dimensional cross-sectional images, and of "Volume Render Mode," a technique to analyze information inside a three-dimensional volume by digitally enhancing individual flow, promises to revolutionize the diagnosis of obstructive sleep apnea patient of level severity of AHI. Using the different post-processing display parameters, the volume-rendered image provides better visualization performance when there are no significant differences in the signal levels of the conventional approach, as shown in Figure   10. In Figure 10, the TKE occurs after critical area for both breathing conditions, either inspiration or expiration, respectively.

Conclusion
In this study, the feasibility of using computational analysis to analyze airflow characteristics has been investigated. Based on the results of this investigation, the following conclusions can be made.