**1. Introduction**

Obstructive sleep apnea syndrome (OSAS) is a potentially life-threatening illness [1–4]. Obstructive sleep apnea (OSA) is a traditional chronic syndrome implicating the adult population, with the highest occurrence reported among middle-aged men [5]. The ailment is characterized by repetitive episodes of a complete or incomplete collapse of the upper airway during sleep, with a consequent decrease of the airflow [6]. Over the decade, several OSA managemen<sup>t</sup> methods have been developed [7,8], and among them by utilizing positive airway pressure (PAP) treatment [4,9], an oral appliance [10] and several nonreversible surgical methods such as mandibular advancement surgery (MAS) [11,12]. MMA is a surgical treatment that involves cutting the upper and lower jaws to realign them [13].

**Citation:** Abdul Latif, M.F.; Ghazali, N.N.N.; Abdullah, M.F.; Ibrahim, N.B.; Razi, R.M.; Badruddin, I.A.; Kamangar, S.; Hussien, M.; Ahammad, N.A.; Khan, A. Modelling the Upper Airways of Mandibular Advancement Surgery: A Systematic Review. *Mathematics* **2023**, *11*, 219. https://doi.org/10.3390/ math11010219

Academic Editors: Mauro Malvè and Lihua Wang

Received: 10 October 2022 Revised: 18 December 2022 Accepted: 23 December 2022 Published: 1 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The improvement of the jaw structures passively persuades an anterior displacement of the soft palate and the tongue, widening the pharyngeal space [12,14]. However, it is noteworthy that MAS is a highly invasive treatment, often associated with complications and aesthetic change [15–17].

Consequently, the treatment should apply to selected patients when all other approaches and first-level surgery have failed or patients with established craniofacial deformities [18,19]. CFD is widely used to fully understand the flow mechanism inside the human airways [20–23]. However, there are contradictions among the studies of the utilization of CFD for OSA study. Therefore, a systematic review is needed to understand the fundamental requirement of OSA CFD analysis fully.

The formulation of this systematic review began with the following research question: How can we properly model the upper airways (UA) as a prediction approach before MAS? Table 1 presents the keyword search string used for this article—the review's process flow, as depicted in Figure 1.

**Database Keywords** WoS ALL = ((maxillomandibular OR mandibular OR mandible\* OR jaw\* OR maxilla) AND (advance\$ OR improvement OR gain OR elevation) AND ("obstructive sleep apnea" OR snoring OR "sleep\* disorder breath\$" OR "pharyngeal airway resist\$" OR "sleep apnea"))

**Table 1.** Keywords for the search string.

**Figure 1.** Flow chart of the review screening process.

### **2. Review Outcome Introduction**

This paper reviews previous studies of OSAS from a CFD perspective focusing on MMA treatment. This review will detail previous CFD technology associated with OSA cases. In this paper, the author will be detailing the rigorous needs of the image processing technique, which is the essential pre-process of CFD modelling.

Furthermore, the UA boundary condition setup will be discussed in detail, starting from the meshing technique of the UA model. The discussion will detail the turbulence model used in the OSA literature, including the setup parameter of UA CFD modelling. Towards the end, the author will discuss related experimental validating processes to indicate the competency of CFD analysis as a tool for MMA treatment assessment.

### *2.1. Airways Imaging Technique*

OSAS has become an interesting topic in research recently. Recent developments in OSA treatment procedures have heightened the need for advanced image techniques [24,25]. The imaging technique available in modern technology has two types: computerized axial tomography (CT) and magnetic resonance imaging (MRI).

One crucial theoretical topic that has engaged researchers for many years is how the UA trigger snoring [22]. The introduction image technique based on acoustics has opened a broad interest in the OSA prediction technique [26]. The MRI or CT approach is confined to a 2D stacking image of the subject's cross-section region, subjected to visualization, essential length, and volume measurement, as shown in Figure 2 [27–29]. The CT or MRI technique, however, has been widely used to evaluate the severity of OSA for surgery determination from the early 90s until today [30–32].

Until today, most of the clinical approaches to OSA treatment depended on the UA's generated cross-sectional image, which provides insight into the UA's narrowing gap and volume [33,34]. The medical decision is solely based on the judgment of the UA crosssection due to time limitation and virtual modeling capabilities limitation [35–37]. These results in non-responsive post-treatment occurred because not all OSA patients responded to surgical treatment positively [38,39]. A more detailed evaluation is needed to properly understand the behavior of OSA in the UA for a more accurate prediction of pretreatment.

**Figure 2.** Example of CT-scan sagittal view of pretreatment of OSA. Adapted with permission from Ref. [28]. 2015, Butterfield.

### *2.2. Modelling of Human UA*

Previously, studies of OSAS depended on statistical reviews of MRI and CT images and various clinical trials [40,41]. The introduction of CFD simulation has revolutionized OSA studies [42]. The knowledge of fluid dynamics is applied to understand the mechanism of OSA in the last decade [43]. The introduction of CFD has allowed studies of modelled UA to explain OSA theoretically [44]. The UA model was first studied based on the teaching

model (Model C12, Carolina Biological Supply Company, Burlington, NC, USA) for medical school students [23] (Figure 3). This model assumes the human airways to be symmetrical. The study by Ted B. Martonen [23] scanned the silicon model into the computer as a 3D model. He then performed the simulation using different flow rates, demonstrating the possibility of applying CFD to model UA.

**Figure 3. Example of** human UA medical school teaching model used by Ted B. Martonen. Adapted with permission from Ref. [23]. 2002, Ted B. Martonen.

In 2003, Heenan et al. [45] developed a 3D model with realistic UA anatomy geometry (Figure 4). This model is an adaptation of the Weibel A model. Their study introduced a CFD method to study the airflow of the UA based on the model, which was less complicated than the actual human UA. Since then, CFD has been an increasingly important area in OSA study.

**Figure 4.** Idealize Weibel A human UA model. Reprinted with permission from Ref [45]. 2003, A. F. Heenan et al.

More literature has emerged that offers contradictory findings of the Weibel A model. Collins et al. [46] compared a geometrically accurate model sourced from the MRI model with the Weibel A model. The result shows a dissimilar comparison flow pattern because the Weibel A model has sharp edge geometry. In contrast, the accurate model has a smooth geometry, as shown in Figure 5. Thus, the simplified model is insufficient to accurately predict human UA's flow behavior. Although the result represented by Collins et al. [46] shows a significant variation in the flow pattern, the model developed by Heenan et al. [45] is a helpful reference for the CFD modelling of UA.

**Figure 5.** Comparison of the Idealize model with the MRI model. Reprinted with permission from Ref. [46]. 2007, T. P. Collins et al.

MRI and CT scan imaging techniques have revolutionized the human respiratory system [47]. This technique enables the construction of accurate geometry of human UA with the help of image processing software such as 3DVIEWNIX, AMIRA, MIMICS, and 3-MATICS [48]. Recent advances in the research of UA flow have emphasized the importance of anatomically accurate UA models. A significant volume of published studies describes the role of MRI or CT in describing the UA flow [49,50]. The first serious discussions and analyses of anatomically accurate UA models emerged during the early 2000s. Xu et al. [51] modelled three cases of OSAS in children aged three to five. The MRI slice image of the children UA was transformed into a 3D model using 3DVIEWNIX software. The study highlights a manual mask filtration technique applied to extract the UA, which removes some of the detail and small voids to simplify the geometry. The simple, clean geometry is necessary to have a good quality CAD model for CFD [52,53].

Most researchers report removing some details of voids and surface smoothing in the UA. It adds complexity to the CAD model, limiting the excellent meshing capabilities and weak CFD convergence [42,53]. Several studies that used MRI or CT show the potential of a functional imaging source for the 3D construction of the UA [25]. However, there are no published data on their sensitivity or specificity [41].

The excellent quality of CAD is a necessity for CFD application. However, the MRI or CT imaging technique requires tremendous effort to construct a good-quality CAD model. A specific technique for creating the UA geometry correct model has not ye<sup>t</sup> been published. The technique is essential because it contributed to the CFD application's dependability and precision in comprehending OSAS [54].

### *2.3. Exclusion of the Nasal Cavity*

Excluding the oral and nasal cavities simplifies the UA model in a significant portion of the UA research. This solved the issue of indefinitely characterizing the wall boundary of the UA. The oral and nasal cavities have less impact on UA flow. Under no circumstances does there seem to be evidence that the nasal cavity may fundamentally alter the pharyngeal flow characteristics, such as the slightest pressure or maximum velocity. The investigation by Zhao et al. [55] demonstrates that a missing nasal cavity will not radically fluctuate the UA's pressure drop and stream velocity profiles. This finding aligns with the examination by Persak et al. [56] and Shah et al. [57]. Calmet et al. [58] conducted a complete study of airways from the trachea up to the start of the nasal cavity (vestibule). In his study, he shows turbulence behavior in the middle of the nasal cavity. This is understandable due to the complexity of the nasal cavity funnel. It is hard to justify the influence of the nasal cavity on UA air flow purely based on a single data sample without comparison with other data in that investigation. Cheng et al. [59] also demonstrated full airways with a nasal cavity comparing pre-surgery and post-surgery data. In his reporting, he did not describe or demonstrate any changes within the nasal cavity, and the ensuing image on the pressure

contour plot demonstrated no substantial change prior to entering the UA funnel. The same study by Ito et al. [10] also showed no significant change in the airflow. Cheng et al. [59] exhibited a full airways model study, including the nasal cavity. However, in that study, the nasal cavity was superimposed. Both pre-surgery and post-surgery provided similar inflow patterns. An investigation by other researchers also revealed the same result [60–62]. The finding clearly shows that excluding the nasal cavity is necessary, as it adds more complexity to the investigation.

#### *2.4. The Meshing of UA Geometry Accurate Model*

As stressed, good-quality geometry is vital for CFD application. A CFD prepossessing step meshing discretizes the geometry into a more distinctive element to make numerical calculation possible. Currently, there are several options for mesh generation available commercially. Typically, as documented in most publications on UA models, most of them use the built-in mesh-generating capability of the widely known ANSYS CFD software [63–66].

Due to the complexity of the geometry, researchers tend to utilize auto mesh generation for retaining the original geometry accurate model of the UA [20,67], resulting in the mesh density concentration on the tight curve or small surfaces. Researchers applied more exceptional mesh cells to solve the unevenly distributed mesh density issue, typically for a UA model consisting of between 500 thousand to 1.6 million cells [55,68]. Typically, for internal fluid flow assessments such as the UA model, a hybrid mesh with at least five inflation layers is used, since it represents the near-wall effect better [20]. Figure 6 shows the example of a hybrid meshing of the UA model.

A finer mesh does not necessarily result in solution convergence. Although convergence is possible with more excellent meshing, it always comes with the high cost of computing. It is undeniable that solution convergence relies on good quality meshing. Most of the UA model research does not report the mesh quality of their model. They depend on finer mesh due to geometry complexity for convergence, which costs unnecessary computing [43]. Researchers utilize the application of the smoothing technique in the hope of achieving better mesh quality [42]. Some of the studies of the UA model utilize third-party meshing software, such as MESHLAB [69], GAMBIT [70,71], and DEP MESH-WORK [72]. Still, most of them only report basic auto mesh generation applications [42]. Recent advancements in the simulation of UA models have increased the demand for improved meshing approaches that result in a higher mesh quality for UA models.

**Figure 6.** Hybrid meshing of a UA model Reprinted with permission from Ref. [43]. 2013, Moyin Zhao et al.

Finer meshing increases the accuracy of the result, which theoretically makes sense as it is close to the source model. However, finer meshing is computationally expensive [73]. In CFD, there is a well-adapted method to minimize unnecessary computing costs. The technique is grid sensitivity analysis [74]. The process determines the optimum mesh density while retaining the accuracy of the result. The study model simulates various element sizes, starting from a coarse to a finer mesh [75]. This method has been well

adopted in the UA model by multiple researchers. Zhao et al. [55] show a plot of axial velocity along a vertical line from the inlet to the larynx wall, starting with a 200k mesh to a most excellent 1.8 million mesh. The result shows that 1.3 million mesh has a similar velocity profile as the 1.8 million mesh while saving 30% of the computing time.

Rahimi-Gorji et al. [73] demonstrate the execution of distinctive grid sizes comprising around 2.4, 3.4, 4.2, and 5.1 million cells and acquire pivotal velocity profiles at two cross areas. They have found that the expansion of grid size measure from 4.2 to 5.1 million cells does not modify the outcomes as it shows an agreeable velocity profile. Subsequently, the study uses a grid with a 4.2 million mesh size for the simulation. However, the author feels that comparing the pressure or velocity profile against the number of cells does not justify the acceptability of the boundary condition. A comparable validation test rig is necessary as a reference to determine the appropriate cell number for the modelling, as shown by Amatoury et al. [76]. A grid size agreeable with accurate experimental data is the ideal mesh size for modelling the UA.

### *2.5. Boundary Conditions*

More literature has emerged that offers contradictory findings of the appropriate boundary condition for UA simulation. The boundary conditions used in UA CFD simulations should characterize and replicate those seen in human respiratory flow. The UA inspiratory flow originates from the nasal cavity and ends at the trachea. In a perfect situation, encompassing static pressure should be characterized at the nostrils and, correspondingly, release airflow at the lower larynx. Numerous published studies define human inspiratory airflow rates differently, as summarized in Table 2. A few examinations determined a time-dependent flow to copy the respiratory cycle.


**Table 2.** Compilation of methods used for CFD modelling of the UA.


**Table 2.** *Cont.*

In contrast, others described a mean airflow stream rate [71,73,77]. Table 2 shows considerable variation in the airflow rate used in UA simulation. The variation is understandable because, ethically, it is difficult to have a human subject measure the actual patient breathing airflow rate. It is also because the human actual breathing flow rate varies from one person to another. Other boundary conditions were customary settings for all the reviewed studies, like a smooth and non-slip wall and a five percent turbulence intensity of the inlet flow [71].

### *2.6. Turbulence Modelling of UA Flow*

The turbulence model adopted in UA modelling varied among these four commonly adapted models: k-epsilon, k-omega, k-omega SST, and large eddy simulation (LES) [63,81]. Turbulent flow features are modelled by solving the variables for kinetic energy (k) and dissipation rate (ε) or specific dissipation rate ( ω)—a literature finding of the adopted turbulence model for UA modelling, as in Table 2. The decision to model a fluid dynamics problem, either laminar or turbulent, is one of the most challenging decisions a fluid dynamicist must take [82]. Therefore, a thorough understanding of the predicted flow field is critical for obtaining the desired outcomes. The estimated Reynolds numbers (Re) range indicates the UA flow to be laminar or transitional. The standard k– ω shear stress transport (k– ω SST) model was proven to be appropriate to simulate this complex flow.

The SST k– ω model has advantages in solving complex transitional flow, including the UA transitional flow [43,83]. The employment of the k– ω shear stress transport (SST) turbulence model provides an enhanced description of flows involving adverse pressure gradients and curved boundary layers [79].

The k–ε model with enhanced wall treatment is suitable for monitoring flow separation with a strong pressure gradient and flow recirculation. These are suitable to determine the flow parameters near the wall of running airflow [63].

Shown by literature, LES has become the preeminent turbulence model for the UA model simulation, as most of the studies show good agreemen<sup>t</sup> with the experimental result [81]. However, a fundamental disadvantage of this turbulence model application is that it demands a high computing expense incurred due to the length of time necessary to solve it. As a result, LES is unsuitable for time-demanding clinical applications. The K– ω SST model also agreed well with the experimental result and is almost comparable with the LES model [43,80]. This turbulence model offers much cheaper computing costs and is suitable for clinical application. Again, the author stressed that an excellent comparison with the physical model provides a justification of the ideal turbulence model for the UA study.

### *2.7. Location of the Inlet*

The velocity profile of many studies demonstrates a brisk airstream at the smallest cross-sectional area at the velopharynx. This flow feature is identified as the 'pharyngeal jet' [46,78,81,83]. Hence, it may be coincidental that the degree of structural difference in UA anatomy systems is virtually tremendous, resulting in a diversity of UA flow patterns.

Before surgery, the highest increase in airflow velocity is observed at the pharyngeal airway when inhaling, and the biggest decrease in pressure is observed downstream, according to Liu's study. Negative pressure and airflow velocity in the whole airway equalized postoperatively. They discovered that velocity reduction at the most constricted portion of the pharyngeal airway correlates most strongly with surgical outcome. However, they emphasized that measurements of non-theoretical dynamic airflow during sleep would be excellent for further validating CFD models [62].

The study by Zhao et al. [55] of flow profiles for both inhaling and exhaling situations indicate that inhaling results in a 30% greater pressure loss and a higher flow velocity than exhaling. The airflow enters the UA stream with higher turbulent kinetic energy and lower total pressure, which expands the likelihood of UA collapse. In short, narrowing the UA increases the velocity of the territorial flow and decreases the pressure. This low-pressure peak may be associated with the severity of OSA, as a high-pressure gradient across the wall border would cause the UA structure to collapse [46,78,81,83]. Furthermore, 10 L/min, 15 L/min, and 30 L/min are the most prevalent input flowrate settings, according to our review [45,73,77–79,81,82].

### *2.8. Experimental Validation*

Experimental validation in the UA model based on recent studies focused on verification checking for numerical studies [43,45,55,64,83,84]. Rapid prototyping is the most common method to create an identically exact UA model. [85]. Heenan et al. [45] have constructed an idealized Weibel A model representing the human oropharynx. The design slightly differs from the one described by Stapleton et al. [86]. The model is two times larger than the original, which describes double the adequate accuracy in setting up using the particle image velocity (PIV) method described in Figure 7. The flow rate of the double scale model was doubled to maintain the full-scale Reynolds number. The construction of the physical model uses the fused deposition modelling (FDM) method with surface smoothing using dichloromethane and an epoxy coating. The result shows a discrepancy between the experimental and CFD. Heenan et al., assume the differences is because of the error in the CFD simulation due to the weak boundary layer of the CFD model [45].

**Figure 7.** PIV experimental setup Reprinted with permission from Ref. [45]. 2003, A. F. Heenan et al.

Collins, Tabor and Young [46] argue that the flow character in the idealized geometry is affected by the shape of the geometries used. The absence of air curvature and surface irregularities in the Weibel A model makes it unusable for predicting the flow pattern. Thus Collins, Tabor and Young [46] developed a hybrid geometry based on the one described by Stapleton, Guentsch, Hoskinson and Finlay [86] and Heenan, Pollard and Finlay [45], as shown in Figure 8. The result are compared with an MRI based model of a 21 years old male. Both were simulated using the CFD method.

**Figure 8.** Hybrid idealized geometry. Reprinted with permission from Ref. [46] 2007, T. P. Collins et al.

The idealized hybrid model shows good agreemen<sup>t</sup> with the one represented by Heenan, Pollard and Finlay [45], as it is a similar geometry. However, the idealized geometry offers a slight advantage over the MRI geometry due to the recirculation zone created by the sharp edge or steps in the idealized geometry. Thus, it is concluded that the idealized geometry is inadequate at predicting the features of the flow of human UA.

Xu, et al. [87] cast a 3⁄4 scale of a geometrically accurate model of UA using a silicon cast for result validation (Figure 9). The reason for using a scaled model is unknown. The study placed eight pressure sensors along the surface of the silicon UA model. The study also placed the same pressure point positions in the CFD model. However, the CFD model did not scale to the cast model. The results agreed to within 2.0 Pa in both inspiration and expiration phases at almost all locations. However, the maximum pressure drop that CFD predicted was 20 Pa higher than in the experiment. The scale difference might play an essential role in pressure distribution; thus, a comparable 1:1 scale model is desirable.

**Figure 9.** 85% scaled silicon cast UA model. Reprinted with permission from Ref. [51]. 2006, Chun Xu et al.

Latter, Zhao, Barber, Cistulli, Sutherland and Rosengarten [43] developed a 1:1 scale model of a geometry-accurate physical model. The model uses the rapid prototyping (RP) method using a transparent polymer. A comparison of pressure distribution along the surface of the UA model is shown in Figure 10. The finding shows good agreemen<sup>t</sup> between the CFD and the experimental result. This finding is significant as it can be used as a base model to verify one's UA model CFD result by following the exact boundary condition and turbulence model.

**Figure 10.** 1:1 scale of the UA model pressure distribution comparison between CFD and the experimental results. Reprinted with permission from Ref. [43] 2013, Moyin Zhao et al.

Due to the complexity and individuality of the geometrically accurate MRI or CT scan model, it is not easy to adjust the geometry to understand the geometrical influence on UA flow better. An idealized geometry should have a better shape agreemen<sup>t</sup> than the average-person UA and not generalize as the Weibel model. The model should have detailed geometry that another researcher efficiently replicates.

### *2.9. Numerical Modelling Issues of UA*

A previous study by Collins has indicated the disadvantage and limitations of the Weibel model. In term of shape based on his study, it is impossible to mimic the organic shape of UA by mean of 3D CAD construction (Idealize model) [46]. The model is oversimplified and has a sharp edge that is not present in the actual human UA. In addition, the model's nasal cavity does not accurately represent the upper airways since it comprises a single large chamber rather than multiple layers of the air passage. The flow mechanism of such a design differs significantly from or completely contradicts the actual UA and specifies the flow characteristic entering the trachea, influencing the UA study's total flow [46]. The authors sugges<sup>t</sup> that the 3D model is the most appropriate geometrical model for comprehending the fundamental flow behavior within the UA. In terms of CAD and meshing, the 3D model is simple to configure for numerical simulation due to its simple shape [45]. CT scans offer an actual or geometry-accurate model of the UA, which can be converted into a CAD model using specialized software [59,63,81]. In numerical modelling, this precise model provides real-world information regarding air movement within the UA. Due to the organic nature of this model, preparing the CAD for CFD analysis necessitated extensive attention to detail [52,53].

Using CT scans to represent UA has advanced the study of UA. However, it is not without flaws. The geometry accuracy can vary depending on the image quality or resolution of the CT scan. This image is difficult to comprehend. Most researchers employ MIMICS to produce the 3D mask of the UA. The masking process requires additional treatment, detailing, or processing, such as filling the unneeded void or removing the surface spike. Poorly processed masks could result in poor meshing and compromise the CFD output [70,71].

The laminar model is unsuitable for UA research because, according to Barber's research, the UA's Reynolds number (Re) ranges from 400 to 3000, indicating that the flow is either laminar or transitional [45,73]. The flow turbulence was modelled using the standard shear stress transfer (SST) k–ω model. Their early sensitivity research applying multiple turbulence models to the physical experimental results demonstrated the suitability of the SST k–ω model in solving the complex UA flow [43]. Therefore, only a few studies have reported the UA utilizing laminar flow [65,71,88]. Compared to experimental results, LES provides the most accurate model of UA due to its superior observation at the wall and

in the center. However, LES simulations require powerful computing capabilities. Other turbulence models, such as k–, are utilized extensively for external flow but infrequently for inside flow [78]. The researcher's findings also indicate that k– has low precision for the UA investigation [88]. The most acceptable and cost-effective turbulence model is k– ω, which offers comparable accuracy to the LES model [81].

Meshing of UA requires hybrid meshing for internal flow so that flow mechanisms close to the wall can be captured accurately. This hybrid meshing requires a high-quality CAD without elements that degrade mesh quality. For the solution to converge, highquality mesh is required [73].
