Simulation of Pharyngeal Airway Interaction with Air Flow Using Low-Re Turbulence Model

Abstract

This paper aims to simulate the interaction between a simplified tongue replica with expiratory air flow considering the flow in the pharyngeal airway to be turbulent. A three-dimensional model with a low-Re SST turbulence model is adopted. An Arbitrary Eulerian-Lagrangian description for the fluid governing equation is coupled with the Lagrangian structural solver via a partitioned approach, allowing deformation of the fluid domain to be captured. Both the three-dimensional flow features and collapsibility of the tongue are presented. In addition, examining initial constriction height ranging from 0.8 mm to 11.0 mm and tongue replica modulus from 1.25 MPa to 2.25 MPa, the influence of both of these parameters on the flow rate and collapsibility of the tongue is also investigated and discussed. Numerical simulations confirm expected predisposition of apneic patients with narrower airway opening to flow obstruction and suggest much severe tongue collapsibility if the pharyngeal flow regime is turbulent compared to laminar.

1. Introduction

Collapse of the human pharyngeal airways during sleep has serious health complications with an estimated 10% of snorers being at risk to obstructive sleep apnoea [1, 2]. These episodes of partial or full cessation of breathing per hour, influence the quality of sleep, reduce brain oxygen saturation and have been linked to hypertension and heart failures [1]. The physiological mechanisms of such conditions are very closely related to flow in compliant tubes or channels. Extensive experiments of flow in Starling resistors have revealed rich variety of dynamics involved in such compressed collapsible tube systems [3, 4].

Not surprisingly, many numerical models have been developed to simulate the system and theorise the mechanisms involved in the collapse and self-excited oscillations of these compliant walls seen in the experiments. Earlier models include lumped parameter models [5] which integrated the flow variables along the whole vessel and various one-dimensional models (e.g., [6–8]) taking into account variation of flow variables in the longitudinal direction along the vessel but assume those variables do not vary within each cross-section. More complicated two-dimensional models have also been developed, such as those by Luo and Pedley [6, 7, 9, 10], revealing certain conditions (function of Reynolds Number and membrane tension), where steady solution are stable and below which self-excited oscillation and vorticity waves are generated in the flow. More recent three-dimensional modelling of collapsible tubes captured strong buckling of the tube [11] and also revealed both possible mechanism for flow-wall energy extraction and critical Reynolds Number for onset of self-excited oscillations [12].

Applying these developments of flow in collapsible tubes to actual physiological phenomena presents a real challenge [1]. On the subject of obstructive sleep apnoea, Chouly et al. [2, 13] and Van Hirtum et al. [14] developed an asymptotic Navier-Stokes equation coupled with linear elastic shell elements to model flow-induced deformation of a simplified tongue subjected to expiratory flow. The collapse of the tongue onto the pharyngeal wall was also simulated via in-vitro experiments to validate their model assumptions and numerical results. In general, narrowing of the airway promotes increased transmural pressure via a venturi effect, resulting in partial collapse of the airway and a nonlinear flow rate retardation as the intraluminal pressure difference is increased—a typical observation in collapsed channel called flow rate limitation.

The coupling between a two-dimensional, quasisteady, laminar, incompressible fluid flow with a shell model developed by Chouly et al. [2] was successful in obtaining efficient real-time results and was validated using pressurised latex in their experiments. Indeed, based on Reynolds Number, the flow is expected to be within the laminar (Re<2000) or transitional (2000<Re<10000) regime. However, the complex morphology of the pharyngeal airway, could give rise to three-dimensional flow features which triggers transition or turbulence regimes at lower Re [15]. In fact, observation by Hahn [16] suggests turbulence kinetic energy at approximately 5% of the inlet velocity is not unexpected in the nasal cavity.

Both neurological and physiological factors have been implicated with apneic syndrome. Previous studies (e.g., [17]) have shown that apneic patients have in general, a narrower opening in the oropharynx region (perhaps, due to obesity or tissue buildup). Pharyngeal compliance among apneic and nonapneic subjects have also been shown to vary significantly, especially during expiratory flow [17].

Hence, the present study intends to instead consider a three-dimensional model of the expiratory fluid flow using a low-Re turbulent model coupled to similar replication of the tongue. The effect of initial airway opening and the tongue stiffness on the collapsibility and flow inside the pharyngeal airway is investigated parametrically. Three-dimensional flow features and turbulence influence are also discussed.

2. Mathematical Model

2.1. Governing Equation for Fluid

Following the approach of Chouly and coworkers [2], a simplified three-dimensional flow configuration representing the pharyngeal airway is illustrated in Figure 1.

Figure 1 

Pharyngeal airway model.

The converging flow through the narrowest opening in the oropharynx (at the base of the tongue) is expected to produce a jet stream and flow separation downstream of the constriction which is characterised by shear instability-induced turbulence [18]. Avoiding the intensive computational efforts involved with a three-dimensional Large Eddy Simulation (LES) and Direct Numerical Simulation (DNS), a more practical approach is employed. Steady, incompressible, turbulent flow in the three-dimensional tube is described using the Reynolds-Averaged Navier-Stokes (RANS) equations coupled with a Shear Stress Transport (SST) 𝑘-𝜔 turbulent model:𝜕𝑢_𝑖𝜕𝑥_𝑖𝜕𝑢=0,(1)_𝑗𝑢_𝑖𝜕𝑥_𝑗=−1𝜌𝜕𝜕𝑥𝑝_𝑖+𝜕𝜕𝑥_𝑗𝜈𝜕𝜕𝑥𝑢^𝑖_𝑗+𝜕𝑢^𝑗𝜕𝑥_𝑖+𝜕𝜕𝜏𝑥^𝑅𝑖𝑗_𝑗,(2) where 𝑢_𝑖 and 𝑝 are the mean- or time-averaged velocities and pressure (𝑖=1,2,3 for three-dimensional analysis), 𝜌 is the fluid density and 𝜏^𝑅𝑖𝑗 is the turbulent Reynolds stress tensor 𝜏^𝑅𝑖𝑗=−𝑢^𝑖𝑢^𝑗,(3) where 𝑢^𝑖 represents the turbulent fluctuating velocity. Using Boussinesq’s assumption, this Reynolds stress tensor is related to the mean velocities by the turbulent eddy viscosity: 𝜏^𝑅𝑖𝑗=−𝑢^𝑖𝑢^𝑗=𝜈^𝑇𝜕𝜕𝑥𝑢^𝑖_𝑗+𝜕𝑢^𝑗𝜕𝑥_𝑖−23𝜕𝜕𝑥𝑢^𝑘_𝑘𝛿_𝑖𝑗,(4) where 𝜈^𝑇=𝑘/𝜔 is the turbulent eddy viscosity.

Therefore, in order to close the RANS equations (1) and (2), only the turbulent viscosity is required, which is determined by using 2 unknowns, that is, turbulent kinetic energy 𝑘 and turbulent rate of dissipation 𝜔. The principle behind SST 𝑘-𝜔 model is to blend the original 𝑘-𝜔 and 𝑘-𝜀 together such that, near the wall, the accuracy and robustness of the 𝑘-𝜔 model is captured, while away from the wall (i.e., free stream or free shear regions), the accuracy of the 𝑘-𝜀 is captured. Thus, the additional two equations required to solve for scalar variables 𝑘 and 𝜔 are [19]𝜕𝑢_𝑗𝑘𝜕𝑥_𝑗=𝜕𝜕𝑥_𝑗𝜈+𝜈𝜎^𝑇_𝑘𝜕𝜕𝑥𝑘_𝑗+𝑃^𝑘𝜕𝑢−0.09𝑘𝜔,(5)_𝑗𝜔𝜕𝑥_𝑗=𝜕𝜕𝑥_𝑗𝜈+𝜎𝜈^𝑇_𝜔𝜕𝜕𝑥𝜔_𝑗+21−𝐹¹𝜎1_𝜔2𝜔𝜕𝜕𝑥𝑘_𝑖𝜕𝜔𝜕𝑥_𝑖+𝛼𝜔𝑘𝑃^𝑘−𝛽𝜔²,(6) where 𝜎_𝑘, 𝜎_𝜔, 𝜎_𝜔2, 𝛼, and 𝛽 are model constants (blended from the original model 𝑘-𝜔 and 𝑘-𝜀 constants), 𝐹₁ represents the blending function which is dependant on the wall distance, 𝑃_𝑘 is the production of turbulence due to viscous forces (function of the mean velocities), and, thus, the last two terms in both equations are effectively representing the production and destruction of each 𝑘 and 𝜔.

In addition, the SST model is further refined by employing a blending of the turbulent eddy viscosity formulation, allowing more accurate prediction of the eddy viscosity near and away from the walls: 𝜈^𝑇=𝑎¹𝑘max𝑎₁𝜔,𝑆𝐹₂,(7) where 𝑎₁ is a model constant, 𝐹₂ is another blending function (similarly depending on the wall distance) and 𝑆 represents an invariant measure of the strain rate [19].

Equations (1), (2), (5), (6), and (7) represent the complete fluid model which is solved by the fluid solver. A low-Re simulation with appropriate near-wall behaviour are readily captured in the SST 𝑘-𝜔  model.

Assuming the flow is symmetrical about its mid-vertical plane, a half-model of the problem is constructed as depicted in Figure 2. Following similar notations by Luo and Pedley [6], the conditions imposed at each boundary are summarised as follows:Inlet: 𝑝=𝑝_in, 𝑥=−(𝐿_𝑢+𝐿/2),Outlet: 𝑝=𝑝_out, 𝑥=(𝐿/2+𝐿_𝑑),Rigid Side Walls: 𝑢_𝑖=0, 𝑦²+𝑧²=𝑑,Rigid Bottom Wall: 𝑢_𝑖=0, 𝑦=−ℎ,Elastic Wall: 𝑢_𝑖=𝑈_𝑖, on Γ, Symmetrical Wall: 𝑢₃=0, 𝑧=0.

(a)

(b)

(c)

Figure 2 

Model definition (a) isometric view, (b) side view looking at symmetrical plane, (c) frontal view looking from the inlet.

No slip boundary conditions at the walls implies that flow velocities 𝑢_𝑖 are zero at the rigid walls and at the elastic segment, 𝑢_𝑖 must match the elastic wall velocities 𝑈_𝑖. At the inlet, a uniform upstream pressure 𝑝_in is specified, which implies a parabolic velocity profile definition at the inlet. While at the outlet, downstream pressure 𝑝_out has been fixed to zero implying a typical stress-free condition.

2.2. Governing Equation for Structure

For a continuum undergoing steady deformation, Cauchy’s equation [20] is effectively a static equilibrium equation: ∇⋅𝜎_𝑖𝑗+𝐟=0,(8) where 𝜎_𝑖𝑗 is the stress tensor and 𝐟 is the external forcing term (𝑖,𝑗=1,2,3 for three-dimensional structures) which in this case, includes external pressure 𝑝_𝑒 plus the fluid pressure and fluid shear at the interface Γ. Mechanical properties of the tongue is not homogenous and anisotropic, more so with varying degrees of muscle activation. For simplicity, a homogenous isotropic material is assumed for the tongue, thus, 𝝈=𝐃𝜺 where 𝜺 is the strain vector and 𝐃 is the constitutive elastic stiffness matrix (which is a function of Young’s modulus, 𝐸, and Poisson ratio, 𝜈 of the structural material).

Considering small deformation theory, second-order terms in the Green-Lagrange strain tensor [20] could be neglected. Thus, the strain-displacement relationship could be described as: 𝜀_𝑖𝑗=12𝜕𝜕𝑥𝑑^𝑖_𝑗+𝜕𝑑^𝑗𝜕𝑥_𝑖,(9) where 𝑑_𝑖 is the displacement tensor. Effectively, (8) and (9) represent the solid model, capturing the elastic wall deformation under external loads.

The tongue anatomy mainly consists of water (with reported density of 1040 kg/m³ [21]) and hence, is generally accepted as incompressible [22]. Accounting for that incompressibility, a Poisson ratio 𝜈 of 0.499 has been used [2, 22] and is also adopted here. Tissues along pharyngeal airway have been reported with a range of moduli: 10–30 kPa for the vocal folds [23], 12–25 kPa for the soft palate [24] and 6 kPa [25] to 15 kPa [22] have been estimated for the tongue. In order to replicate the response of a real tongue, a hydrostatically pressurized latex shell tube was utilized by Chouly et al. [2, 13]. Similarly, a simulated shell modulus 𝐸 in combination with an external pressure 𝑝_𝑒 that mimics this response is followed in this paper.

The elastic wall is assumed fixed everywhere except at the common face interfacing with the elastic segment of the fluid domain where an external pressure 𝑝_𝑒 is imposed simultaneously with the pressure 𝑝 applied from the fluid domain.

3. Computational Method

A commercial finite volume solver (CFX) is used to solve the RANS equations in (1), (2), (5), (6), and (7). The fluid domain is subdivided into 8-noded three-dimensional hexahedral elements leading to an assembly of discretised algebraic equations in terms of the nodal unknowns 𝑢_𝑖, 𝑝, 𝑘, and 𝜔. In general, finer meshes are located closer to the walls (elastic, rigid side and bottom) where higher velocity gradients are expected and flow separation or reattachment normally occur. In addition, denser meshes are also located in the vicinity of the constriction where sharp changes in pressure are expected.

The governing fluid equations are based on a fixed Eulerian frame of reference in space. In order to account for boundary deformation and thus, deformation of the fluid mesh, an Arbitrary Lagrangian-Eulerian (ALE) description is employed in CFX. Effectively, in an ALE framework, the mesh velocity 𝑢_𝑗̃ is subtracted from the material velocity 𝑢_𝑗 in the differential operator [20, 26] on the left-hand side of transport (2), (5), and (6). Thus, an additional nodal unknown (i.e., mesh velocity 𝑢_𝑗̃) needs to be solved. CFX employs a diffusion model, which smoothes the internal domain nodal displacement according to a Laplace differential equation [19]: ∇⋅Γ_disp∇𝛿=0,(10) where 𝛿 is the nodal displacement relative to previous mesh positions, Γ_disp is some mesh stiffening parameter and (10) is solved subjected to specified nodal movement at the boundaries.

The commercial finite element solver (ANSYS) is used to solve the partial differential equations (8) and (9). Similarly, the structural domain is discretized into 8-noded SOLID185 elements and (9) is approximated using interpolation functions in terms of the unknown nodal deformations. Minimization of the variational total potential energy or weighted Galerkin residuals lead to an assembly of algebraic equation which could be solved as a function of imposed boundary conditions.

The fluid-structure interaction is achieved by satisfying either a velocity or displacement continuity (i.e., (11) or (12), resp.) and force equilibrium (13), at the common interface between both fluid and solid domain:𝑢𝑠_𝑖=𝑈_𝑖onΓ,(11)_𝑖−𝑝+𝜏=_𝑖𝑗𝑑_𝑖onΓ,(12)=𝜎_𝑖𝑗onΓ,(13) where 𝑢_𝑖 and 𝑠_𝑖 are, respectively, the flow velocity and boundary displacement on the interface, 𝑈_𝑖 and 𝑑_𝑖 are the structure velocity and displacement, respectively, on the interface. The terms on the left-hand side of (13) represent the fluid state of stress and the right-hand side is the stress on the structure.

In order to match these conditions simultaneously in both fluid and solid solvers, a successive iteration is employed in an ANSYS Workbench platform. The fluid variables are solved based on the initial geometrical configuration. The fluid pressure and shear forces at the elastic segment are interpolated to the nodes on the elastic wall. Applying them together with the external pressure and imposed boundary conditions, the solid solver solves for the deformation of the elastic wall. These nodal deflections are then interpolated back to nodes on the elastic boundary of the fluid, which is used to effect the mesh deformation on the fluid domain. The fluid solver then solves the unknown fluid variables using the current geometrical configuration and the process is repeated until the nodal deflection and forces in (12) and (13) from current and previous coupling iterations are within a specified tolerance. The overall process can be summarised as in Figure 3.

Figure 3 

Fluid-structure coupling flowchart.

Note that for steady conditions (where inertial effects are neglected and structural deformation reaches static equilibrium), elastic wall velocities 𝑈_𝑖 approaches zero.

Several runs with different downstream length 𝐿_𝑑 suggest that numerical results are not sensitive to the outlet proximity. In addition, finer mesh discretization also indicates similar overall wall shear distribution and more importantly, similar separation location off the elastic wall, suggesting adequate discretization.

4. Results and Discussion

In order to validate the structural model, a Poisson ratio of 0.499 and a Young’s modulus 𝐸=1.75 MPa for the isotropic elastic wall was used, giving similar load-deflection response of the tongue replica to external pressure 𝑝_𝑒, as expected in [2] (Figure 4).

Figure 4 

Pressure versus deflection curve for elastic wall. (Deflection is measured as lowest point on the symmetric plane. Normalization constants 𝑃𝑜=1000 Pa and ℎ𝑜=1.6 mm.)

4.1. Parametric Investigation

For the purpose of studying the influence of the following physiological variations on the collapsibility and flow pattern in the pharyngeal airway, the following range of parameters were simulated:(i)initial constriction height (ℎ)0.8-11.0 mm,(ii)elastic wall modulus (𝐸)1.25-2.25 MPa.

The initial constriction height is varied to simulate the effect of different degrees of opening behind the tongue within the oropharynx in apneic or nonapneic subjects. The opening gap behind the tongue is typically 11 mm [27] but is perhaps lower when subjects are in sleeping position. Preliminary study by Chouly et al. [18] reported an opening of between 1 to 2 mm during sleep supine position. In line with this, variations between 0.8 to 11.0 mm have been investigated in this paper.

The influence of the modulus of the tongue is also of interest. Variation of mechanical properties in the human population is not unexpected. The degree of variation in pharyngeal compliance measured by Brown et al. [17] among apneic patients could offer some insight into the variation of stiffness in the retroglossal region. Variations in the tongue stiffness is also indicative of the degree of muscle activation [21] and have been used to capture the directional activation of muscle tissues. In this paper, considering the replication of the tongue by an elastic shell, a±0.5 MPa variation is considered.

In order to investigate the first parameter, the elastic modulus 𝐸 and external pressure 𝑝_𝑒 was fixed to 1.75 MPa and 200 Pa, respectively, while the initial constriction height ℎ_𝑜 was varied. The nonlinear variation of the flow rate with intraluminal pressure difference Δ𝑝=𝑝_in−𝑝_out is shown in Figure 5(a) for various ℎ_𝑜, revealing typical flow rate limitation phenomena with increased expiratory efforts. In particular, the increased flow rate and reduced degree of flow plateauing, suggest much reduced flow obstruction or collapse with increasing ℎ_𝑜. In addition, the lower degree of flow plateauing for the ℎ_𝑜=11 mm case would suggest a near-rigid flow rate is achieved and demonstrates the relatively low obstruction in nonapneic patients.

(a)

(b)

Figure 5 

Parametric investigation (a) effect of initial airway opening ℎ𝑜 on flow rate 𝑄 and constriction height ℎ for a fixed tongue replica modulus 𝐸=1.75 MPa, (b) effect of tongue modulus 𝐸 on flow rate 𝑄 and constriction height ℎ for a fixed initial opening ℎ𝑜=1.2 mm (𝑄𝑜=10 L/min, 𝑃𝑜=100 Pa, and ℎ𝑜=initialconstrictionheight).

The second parameter investigated is the effect of tongue modulus on the expiratory flow interaction with the tongue. In order to investigate this parameter, a configuration with initial constriction height of 1.2 mm and 𝑝_𝑒=200 Pa was used. The effect of softening and stiffening of the tongue modulus is summarized in Figure 5(b). Reduction in modulus translates to increased compliance of the airway cross-section to the applied pressures (as suggested by reduced constriction height ℎ in Figure 5(b)). Thus, exarcebating its susceptibility to collapse and increases resistance to flow rate. The approximately linear trend of collapse with increasing Δ𝑝 is consistent with the venturi effect, where the reduction in flow area underneath the elastic wall increases the flow velocity, thus, generating much lower mean pressure (suction) that further collapses the wall.

Comparing the relative influence of both parameters, ℎ_𝑜 and 𝐸, to the wall collapse and flow rate, suggest that both parameters are important. For instance, a 40% stiffening of the tongue modulus (𝐸=1.25 to 1.75 MPa) leads to approximately 11% reduction in wall collapse, as opposed to 10% reduction in collapse for a 50% increase in initial airway opening (ℎ_𝑜=0.8 to 1.2 mm).

Based on the hydraulic diameter, 𝐷=4𝐴/𝑃 where 𝐴 is taken as the cross-sectional area at the inlet and 𝑃 is the inlet perimeter, the Reynolds number, Re simulated in this investigation ranges from 274–9520. Considering a transitional flow regime in the pharyngeal airway replica (by using a low-Re turbulent model), similar to those observed in previous laminar investigations, narrower opening along the airway predisposes patients to larger wall collapse and hence, more severe flow rate limitation.

4.2. Elastic Wall Deflection

The cross-section along the narrowest opening (for a ℎ_𝑜=1.2 mm case) is shown in Figure 6 for several intraluminal pressure differences Δ𝑃. As expected, in addition to increased deflection as Δ𝑃 is increased, the profile also reveals a larger area close to the side walls where the tongue replica is supported and tends to narrow nonlinearly towards the symmetric plane 𝑧=0. This profile perhaps influences the flow patterns and will be discussed in the next section. The difference in 𝑦-position at the side wall (𝑧=0.0125) is due to the different downstream location of narrowest opening. This is associated to the downstream deflection of the elastic wall as pressure is applied in the axial direction to the upstream portion of the wall.

Figure 6 

Deflection profile in the lateral direction across narrowest opening (see Figure 2(c)).

In comparison to laminar flow, the low-Re turbulent flow predicts a more severe elastic wall collapse, as illustrated in Figure 7. This would suggest a much severe pressure drop is experienced by the elastic wall in transitional or turbulent flow as oppose to laminar flow. Indeed, examining the pressure contour in Figure 8, minimum suction pressure of −81 Pa underneath the constriction is predicted in the low-Re turbulent flow compared to approximately −30 Pa in the laminar flow. This seems to be consistent with numerical findings on a rigid pharynx by Shome et al. [15] where at similar flow rates, a 40% increase in pressure drop is predicted in turbulent flow compared to if the flow is laminar and this is expected to be more severe if the flow boundaries are collapsing, as indicated here. Thus, flow regime in the pharyngeal airway has significant effect on the degree of apneic obstruction and is important to be considered in pre-operative treatments.

Figure 7 

Elastic wall deflection with increasing Δ𝑃 for a ℎ𝑜=1.2 mm case. (Normalization constants: 𝑃𝑜=100 Pa and ℎ𝑜=initialconstrictionheight.)

(a)

(b)

Figure 8 

Pressure distribution in the vicinity of the constriction at Δ𝑃=200 Pa for case with initial constriction height = 1.2 mm and 𝑝𝑒=200 Pa (𝑃𝑜=100 Pa, 𝜈𝑜=10 m/s and ℎ𝑜=15 mm [2]).

4.3. Three-Dimension Flow Features

Typical flow patterns in the idealised airway are captured in Figure 9. In general, flow separates from the elastic wall, the bottom wall and the side walls. Three-dimensional recirculation is observed immediately downstream of the elastic wall above the jet streams in all Δ𝑃. Not surprisingly, swirling about the 𝑥-axis is also formed as the jet separates off the rigid bottom wall and is forced to follow the circular cross-section of the idealised airway.

(a)

(b)

Figure 9 

Streamline plot and contours of velocity along vertical symmetry plane, bottom wall and several axial cross-sections for (a) Δ𝑃=20 Pa and (b) Δ𝑃=200 Pa (ℎ𝑜=1.2 mm case).

Streamlines from upper and bottom half of the inlet reveals the tendency for the stream to migrate towards the side walls, where the cross-sectional area is more open. This side stream (as is also visibly shown by the velocity contours peaking towards the side) is more pronounced at higher Δ𝑃 where the elastic wall deflection is further skewed downstream. Thus, the swirling strength of the inner recirculation core is increased as the outer streams follow the curvature of the walls.

Consequently, the pressure profile is also influenced by the lateral distribution of the cross-sectional opening and velocity. This is evident from the pressure profile plots along axial planes located at several lateral locations shown in Figure 10 where the critical suction pressure in the vicinity of the constriction is somewhat reduced closer to the side walls where a larger opening is anticipated.

Figure 10 

Pressure profile along several axial planes at top wall.

The influence of the side wall onto the lateral pressure and velocity profile is also expected to influence the location of flow separation off both the tongue replica (elastic wall) and oropharynx base (rigid bottom wall). Figure 11 presents the axial shear profile along the bottom walls at several lateral locations 𝑧=0, 𝑧=0.005 and 𝑧=0.0, revealing a shift in the position of zero wall shear as the flow approaches the side walls. Thus, suggesting a delay in flow separation as the side wall is approached.

Figure 11 

Wall shear profile along several axial planes at the bottom wall (note that elastic wall is located within −0.02136m<𝑥<0.02136m).

5. Conclusion

A three-dimensional low-Re turbulent rendering of the airway-tongue replica investigation by Chouly and coworkers [2, 13] and Van Hirtum et al. [14] have been conducted. The low-Re turbulent fluid-structure interaction simulation predicts a more severe pressure drop in the constriction compared to a laminar fluid-structure interaction, reflecting a similar trend as the turbulent rigid pharynx simulation by Shome et al. [15]. Similar flow rate limitation is observed with laminar flow, although more pronounced narrowing of the area behind the tongue is generated.

Parametric study representing some geometrical and mechanical factors that affects apnea in the airway was investigated. The investigation suggest that both a narrower opening and an increased compliance of the airway increases the susceptibility of apneic episodes. In addition, three-dimensional flow features relevant to the airway-tongue replica is also presented.

Brown et al. [17] had successfully measured changes in cross-sectional area along the pharyngeal airway during expiratory and inspiratory breathing pressures, on both apneic and nonapneic patient using an acoustic reflection technique, allowing pharyngeal compliance to be estimated as a function of the distance from the mouth. For future work, this varying stiffness along the pharyngeal airway could be incorporated in order to capture the interaction between multiple segments of differing compliance within the pharyngeal airway and how it influence collapse and reopening of the lumen. Ultimately, a real pharyngeal airway model could be used for the fluid-structure investigation.

Acknowledgments

The authors acknowledge financial support from the Government of Malaysia via the sponsorship by the Ministry of Higher Education under the IPTA Academic Training Scheme awarded to the first author. The financial support provided by the Australian Research Council (project ID LP0989452) and by RMIT University through an Emerging Researcher grant is also gratefully acknowledged. In addition, many thanks to Dr. Franz Chouly for many fruitful communications.