Biotransport in human phonation: Porous vocal fold tissue and fluid–structure interaction

The human voice is produced at the glottis within the larynx by flow modulation that results from the two oscillating vocal folds (VFs). More precisely, the pressurized air from the lungs creates a pulsatile jet in the larynx and interacts with the VFs to initiate flow-induced vibrations. Therefore, the phonation process involves complex biomechanical fluid–structure interaction (FSI) between the VFs and the glottal airflow in the larynx, where the tissue vibrations and fluid dynamics are coupled. The majority of voice disorders arise due to VF structural abnormalities and have a one-to-one correspondence with these FSI features of phonation.¹ For example, vocal fatigue is caused by the high subglottal pressure during the vibratory pattern, the possible cause for hoarseness is aperiodic vibration, aphonia is caused by an inability to produce VF self-oscillation, and a breathy voice is the result of incomplete glottal closure.² 

Voice disorders have a prevalence of almost 30% in the general population³ and 60% in professions with high voice usage, such as educators, public speakers, and singers.⁴ It has been estimated that between 3% and 9% of the population in the United States suffer from some type of chronic voice dysfunction that significantly affects their career and daily communications.^5,6 Thus, understanding fundamentally the FSI in the speech process is essential for studying VF pathology and the physics of voice, enhancement in the prevention of voice disorders,¹ and developing computer-based tools for their clinical management.⁷ 

Computational modeling of phonation provides a prediction of VF dynamics and glottal flow, which can lead to an improved understanding of their interactions. Current efforts to develop computational models of phonatory biomechanics are critically reviewed in Ref. 8, including pure fluid flow models to study changes in aerodynamics within the larynx while neglecting structural vibration of the tissue^9,10 and FSI studies of VFs.^11–16 However, a major modeling limitation in some of the previous research was inducing the vibration by applying direct load pressures to the VF surfaces to study the structural tissue dynamics,^9,10 which fails to provide insight into the realistic physics of phonation due to neglecting the fully coupled behavior of airflow and tissue. Moreover, while substantial progress has been made in studying FSI models that couple glottal aerodynamics with tissue vibration, considering the folds as an elastic structure is a common assumption in all of them. More precisely, the porous structure of the VF tissue is ignored in the previous FSI simulations available in the literature,^4,14,16,17 while in reality, the VF tissue is a biphasic material.¹⁸ Indeed, inside the VF mucosal covering, the tissue is loosely composed of specialized proteins, carbohydrates, lipids, collagen, and elastin fibers, while the liquid such as extracellular fluid flow fills the spaces among these elastic solid structures.

Essentially, since most soft tissues are fully saturated with water, poroelastic models are suitable to describe their deformations. Porous media equations have been used for modeling the transport of drugs and lipids in blood vessels,^19–22 blood flow in the myocardium,^23,24 intradermal injection through skin tissue,²⁵ as well as the interstitial fluid in articular cartilage,²⁶ and intervertebral disks.²⁷ Poroelasticity has also been applied to analyze the mechanical properties of the brain.^28,29 The behavior of the poroelastic materials has been studied extensively; however, only a few studies have included poroelastic models in FSI simulations, see for example, Refs. 21, 30, and 31 and references therein. In these studies, the biological tissues are modeled as biphasic with the properties governed by the porous solid and the fluid occupying the pores which significantly improve our knowledge of the mechanical properties and microstructures of the tissue. Particularly in the context of human phonation modeling, the 80% fluid volume fraction of the VF tissue justifies the use of poroelastic models.³² Considering a biphasic model of the VF allows for the inclusion of systemic hydration (i.e., the fluid within the VF body) of the tissue,³³ which enables greater accuracy in the computational phonation models since the hydration has a significant impact on the function of the tissue and affects its viscosity and elasticity.³⁴ 

Hydration also has an important role in the prevention and treatment of voice problems, as well as in the ease of voice production. Namely, local changes in tissue perfusion is often considered an important biomechanical factor contributing to vocal fatigue.³⁵ The theory of a porous VF model can also provide additional information on the fluid circulation within the tissue. Some voice disorders resulting from localized VF lesions (such as polyps and nodules³⁶) are related to lower perfusion levels.³⁷ In fact, reports have suggested that blood flow to the VF and tissue perfusion fluctuate during phonation.^38,39 This is in agreement with observations in Ref. 40, which indicated that the mild inflammation after the vocal loading task decreased oxygen pressure in VFs and induced hypoxia. In another experimental study on rabbits, it has been observed that over-vibration of the VFs due to voice overuse or abuse significantly affects the permeability of the tissue.⁴¹ However, despite the scientific importance, the impact of VF permeability on filtration flow and the extent that phonation conditions affect the flow to the tissue has not been studied before. There is, thus, a profound need for a greater understanding of the fluid flow in the permeable VFs and investigating its contribution in phonation models.

A number of published reports have made progress in this direction and included the porous structure of the VF tissue to compare the phonation process according to permeable features of the human VF^35,42–47 where Darcy's law governs the movement of fluid through a porous medium. The poroelastic theory has been used to describe the liquid redistribution in the VFs during phonation,⁴² as well as to explore the relation between hydration, tissue permeability, and vibration amplitude.⁴³ However, the glottal aerodynamics is ignored in all these previous studies that have considered the biphasic model of the VFs, and the vibration was imposed by applying pressure loads to the fold surfaces,^35,44,45 instead of a flow-induced vibration. More precisely, these aforementioned models do not use self-oscillation to drive and sustain phonation; instead, boundary loads in the form of distributed pressure were applied directly on the subglottal and intraglottal surfaces of the fold to create the vibration pattern which neglects the FSI features. Moreover, porous VF tissue has been modeled in Refs. 44 and 45 with simplifications that do not reflect the true geometry of the fold such as a three-dimensional rectangular parallelepiped,⁴⁵ or by considering only onefold based on the symmetry assumption about the glottal midline which does not allow for possible alterations in the symmetry of the problem.⁴⁴ Therefore, previous computational studies of porous VFs fail to simulate the realistic physics of phonation due to neglecting the complexity involved in resolving the fully coupled FSI interaction between a poroelastic tissue and the glottal airflow or by disregarding the incorporation of a realistic VF geometry.

In order to overcome the shortcomings of previous porous VF models, FSI must be employed in a real-life model of larynx to describe various biomechanical tissue characteristics associated with the coupled physics of the phonation. This is making the system a fluid–porous structure interaction (FPSI).¹⁸ While each of these components, namely tissue porosity and FSI, are well-studied using numerical simulations, there are no mathematical models or computational studies that capture both the flow inside a VF and its interaction with airflow dynamics. Incorporating this integration is the white space that the current paper aims to cover. The previously mentioned knowledge gap is because coupling between a fluid and porous structure is a computationally challenging problem, utilizing substantial computational resources. However, the phonation process is a strongly coupled flow–structure interaction problem, and oversimplification of the glottal flow physics may not yield reliable results in determining glottal pressure forces and the resulting VF dynamics during self-oscillation. Incorporating a FPSI approach is necessary to accurately represent the true nature of the VFs and to provide a more realistic model to understand VF behavior. Moreover, although the aforementioned porous studies were found to suggest significant attribution of filtration flow in VFs, no previous research has attempted to explore the influence of lung pressure and VF permeability on it, using an integrated FSI methodology capable of simulating the complex interactions present in this coupled system. The limited attempts to study tissue hydration by coupling glottal airflow with poroelastic VF tissue (such as in Ref. 48) were unsuccessful in coupling the full poroelastic formulation with airflow equations, and therefore, interstitial fluid movement within the tissue was estimated from the hydrostatic stress gradient.

The main objective of this study is to develop the coupled glottal airflow–porous VF structure computational framework to quantify the interstitial fluid dynamics within the tissue and identify the association of fluid circulation with porous properties and lung pressure. To investigate the effect of phonation conditions on interstitial fluid distribution and determine the functional importance of porous structure in this regard, tissue's permeability coefficient and inlet lung pressure are varied and the resulting model predictions, including filtration velocity distributions within the VF tissue are evaluated. We believe that this work is the first which specifically investigates the interaction between the laryngeal aerodynamics and vibrating porous VF tissue.

The manuscript is further organized as follows: Sec. II covers the problem formulation including the geometrical configuration and governing equations for the representation of the glottal flow, VF model, and filtration flow. It also summarizes the computational setup to perform the simulations. The computational results for a baseline case on fluid flow data including the pressure and velocity fields, as well as tissue deformation metrics are covered in Sec. III. Additionally, a parametric study is performed, and the results are compared for simulated VF cases. Finally, we summarize the significant findings from this work, and the conclusions and future directions are presented in Sec. IV.

A sketch of the human larynx and VFs in the coronal plane is displayed in Fig. 1. We study the problem of the interaction between the glottal flow (⁠ $Ω f$⁠) with deformable VFs made of a water-saturated poroelastic material (⁠ $Ω p$⁠). A schematic of the three-dimensional idealized model of the larynx containing VFs under consideration is presented in Fig. 2. The VFs and the glottal airflow regions are shown in blue and gray respectively, while the filtration of fluids through the porous VF media is identified by curved arrows [Fig. 2(a)]. The vibration-induced fluid movements inside the porous tissue as well as the fluid exchange through the lateral boundary simulating water resupply from the blood vessels are considered. The sketch in Fig. 2(b) demonstrates the interconnected porous network inside the tissue, along with an illustration of the lateral, medial, superior, and inferior surfaces. The curved black arrows show the circulation of fluid inside the tissue, while the red arrows represent the fluid exchange between the blood vessel and VFs.

The rectangular channel which represents the air domain is divided into subglottal, supraglottal, and glottal regions [Fig. 2(a)]. The sub- and supraglottal lengths are selected based on previous experimental results⁵¹ to capture the whole glottal flow. The cross-sectional geometry of the VF domain is defined by aligning with the geometric features of a laryngeal CT scan of a human subject in Ref. 52. The shortest distance between the two VFs along the medial surfaces is time-varying during the VF deformation and is defined as the minimum glottal gap width, denoted by $d G$⁠, with the minimum distance of 0.02 cm.⁵³ 

Table I summarizes geometric dimensions of the laryngeal-VF model as displayed in Fig. 2 and the corresponding values of these parameters. The considered geometry is comparable with the widely used simplified M5 model proposed by Thomson et al.⁵⁴ and Scherer et al.,⁵⁵ and also corresponds to the larynx model in Ref. 56. The VF geometric parameters include the length (⁠ $L V F$⁠) in the anterior–posterior (AP) direction, thickness (⁠ $T V F$⁠) in the inferior–superior direction, and depth (⁠ $D V F$⁠) in the medial–lateral direction of the VFs. These morphological dimensions are chosen within the human physiological range measured by previous studies and relevant clinical data.^57–59 In particular, the values for the $L V F$⁠, $T V F$⁠, and $D V F$ in Table I represent the average between two values based on the upper and lower range for male and female morphological data for each parameter. The length of the VF is found to be in the range of 0.9^60,61–1.5 cm,^58,62 while its thickness and depth values vary with ranges of 0.5⁵⁷–1.1⁶¹ and 0.5⁵⁷–1.2 cm,⁶³ respectively.

The flow direction is along the axial midline, and the computational domain is assumed to be symmetric about the glottal channel centerline. The VFs, glottis, and shapes of the false VFs comprise the glottal space. The FSI interface of VFs (see Fig. 2) is subject to the influence of the intraglottal airflow pressure. The subglottal airflow pressure is identified with the pressure at the inlet, $P i n$⁠, and the airflow is driven from the inlet (left) to the outlet (right).

In Sec. II B, we summarize the governing equations to simulate the aerodynamics, VF fluid circulation, and the mechanics of the VF models.

Here, $v$ and $ρ f$ represent the fluid velocity field in the larynx and air density, respectively, and $σ f=− p fI+2 μ fD(v)$ stands for the fluid Cauchy stress tensor where $p f$ is fluid pressure and $μ f$ denotes air dynamic viscosity. The symmetric part of fluid velocity gradient is defined as $D v= 1 2(∇v+∇ v T)$⁠. The airflow in the larynx can be regarded as incompressible Newtonian fluid for low Mach numbers (Ma < 0.3), as is the case in human phonation. Therefore, the fluid is assumed to have the constant density of $ρ f$ = 1.185 kg/m³ and viscosity of $μ f$ = 1.831 × 10⁻⁵ kg/m s, representing airflow properties at 25 °C room temperature. Contingent on the computed Reynolds number in Sec. II C, the flow is considered to be laminar, similar to many previous models, which provides reasonable approximation for laryngeal flow (see Ref. 9 for detailed discussions).

In Eq. (4), $E s$ and $μ s$ represent the Young's modulus and the Poisson's ratio of the tissue, respectively. These stiffness parameters are considered to be $E s$ = 40 kPa and $μ s$ = 0.45, similar to previous studies,^65,66 and experimental measurement.⁶⁷ The symmetric part of the deformation gradient for the tissue is denoted by $D U= 1 2(∇U+∇ U T)$ under the assumption of small deformations. It should be noted that considering the VF as a one-layer isotropic, linear elastic material is a common simplifying assumption in models of airflow and VF tissue interaction.^68,69 Justification for linear analysis has been elaborated in Ref. 70, and it has been claimed that during phonation VFs exhibit a nearly linear stress–strain relationship when active muscular tension is present.⁷¹ In particular, linear elasticity was an approximation because the strains observed in the model fall within the linear region.⁷² Moreover, although histological studies have indicated that VF has a multi-layered structure,⁷³ it has been observed to behave mechanically as a one-layer structure for most phonation conditions.⁷⁴ 

The material properties of the VFs and airflow used in numerical simulation of the baseline case are summarized in Table II. Material parameters are chosen based on those reported in the literature.^75,76 The permeability is set to be $K V F$ = 10⁻¹³ m², which is within the range measured in a bovine VF scaffold.⁷⁷ 

As outlined in Secs. II A and II B, the computational model consists of two domains: the air in the larynx and the VFs. Glottal airflow is considered unsteady and incompressible, VFs are assumed to be biphasic with liquid-filled pores and a homogeneous and isotropic elastic skeleton, and the fluid flow through the connecting pores of the porous medium is modeled using the Brinkman equations. The computational 3D model of the larynx and boundary conditions is shown in Fig. 2.

The following assumptions are further made to simplify the physics involved. First, a one-way interaction between the interstitial fluid and VF structure is considered. Specifically, the fluid motion within the VFs is influenced by VF vibration, but the influence of fluid motion on VF vibration is ignored. More precisely, it is assumed that the elasticity of the structure does not play a significant role on the calculated porous pressure and the corresponding filtration flow. However, it should be noted that two-way interactions (fully coupled FSI) between the glottal airflow and elastic VF deformation is assumed to compute the realistic glottal fluid pressures acting on the VFs. While the glottal flow, instead of the elastic property of the material, is considered the core component of the filtration velocity in the VF model, the simplified problem still retains the main difficulties associated with the fluid-porous media coupling. Second, the first term in Eq. (5), i.e., the time derivative $∂ q p ∂ t$⁠, is extremely small. Hence, it is disregarded. The effect of gravity on both the solid and fluid phases is also neglected. Finally, water exchange from the blood vessel is assumed to occur only across the lateral boundary surface.

Dirichlet boundary conditions are applied for the pressure at the inlet and outlet of the fluid domain, with values consistent with physiologically realistic measurements by Holmberg et al.⁸⁰ A static inlet pressure of $P i n$ = 1 kPa is defined at the subglottal region while the pressure at the outlet of the supraglottal region is fixed at $P out$ = 0 kPa. At the inflow and outflow, a zero normal velocity gradient condition is applied across the domain, which allows the flow rate to settle to a value dictated by the pressure drop and the flow impedance. No-slip boundary condition is defined on all bounding surfaces of the fluid domain except the inlet and outlet. The fluid–solid interface is defined between airflow and the superior, medial, and inferior surfaces of the VFs (see Fig. 2). These FSI surfaces of the VFs in contact with air are allowed to deform and move with no-penetration conditions, whereas the other wall regions are assumed rigid. More precisely, a fixed boundary condition is applied to the VF lateral, anterior, and posterior surfaces that are in contact with the arytenoid cartilage. The lateral surfaces of the porous VFs in contact with the outside are treated with opening conditions to model the fluid exchange between the blood flow and the tissue which allows for fluid to circulate within the folds (see Fig. 1, red arrows). This means that the pressure value, $p p$⁠, on this boundary will be calculated by the solver. It is assumed that the filtration flow is fully developed on the lateral boundary. Moreover, it is assumed that there is no relative movement between liquid and solid in the anterior and posterior boundary surfaces (i.e., no-penetration condition). The airflow is coupled to the VF tissue deformation problem by prescribing two conditions at the common air–tissue interfaces: the continuity of the normal stresses along the FSI surfaces is enforced, and the displacement of the interface for the fluid and VF tissue are prescribed to be equal which indicates that the fluid adheres to the structure. Finally, contact between the two folds is modeled by creating a frictional contact between the medial surfaces of two VFs and enforcing the initial glottal gap of 0.02 cm as the minimum distance between the nodes of the contact pair during the simulation. This value is in the range reported by Ref. 53 and has been used in previous phonation models.^81,82 Imposing this gap instead of a complete glottal closure is necessary to prevent the flow domain from being disconnected during the VF closure and to ensure the success of the flow solver.

The parameters, conditions, and equations listed in Secs. II A–II C are implemented in the commercial code ANSYS® Workbench (version 2022.R2, ANSYS Inc., Canonsburg, PA, USA) software to carry out the simulations. Aerodynamics and interstitial flow equations are modeled and solved using the computational fluid dynamics (CFD) software package in the ANSYS CFX component to capture the velocity and pressure of the air traveling through the larynx as well as the perfusion within the VFs, while the deformations and mechanics of the VF tissue are defined in the ANSYS mechanical component. The CFX and mechanical components are joined and interact via a single FSI boundary, whose behavior is dictated by the ANSYS Workbench platform System Coupling module.⁸³ The System Coupling component solves the fluid and solid domains separately to first obtain pressure loads from the airflow in CFX that are applied to the VFs on the fluid–solid interface, solving for tissue displacement in Mechanical. This process of the System Coupling algorithm in ANSYS is described in detail later in this section.

Furthermore, the fluid circulation within the permeable VF tissue is simulated in a second ANSYS CFX component that receives the glottal fluid flow forces of the FSI boundary as a general momentum source term. More precisely, the FSI loads obtained from solving the coupled glottal fluid flow (CFX) and VF deformation (Mechanical) system are exported from the CFX component in the form of pressure gradient. These data are imported into the second CFX component that contains the permeable VF model and are applied to the VF FSI interface to capture the filtration velocity and pore pressure within the porous tissue. In order to apply the FSI load (i.e., pressure forces of the airflow) on the FSI boundary of the porous VF, a subdomain is created in CFX which includes the whole VFs domain, and a general source term is added to the momentum equation in the aforementioned subdomain. The specified momentum source is set as an expression function of the pressure gradient variable, and an interpolation function within CFX is used to create a function from a set of coordinates paired with the values of the pressure and the resulting pressure gradient variable at those coordinates. More precisely, the user function provides forces in three dimensions on the FSI boundary to capture the vibration-induced filtration flow in the porous VFs. Through this process, the three components of FSI force in Cartesian are transferred to solve for the pore pressure and filtration flow.

The process of the System Coupling algorithm in ANSYS is as follows. The Navier–Stokes and continuity equations governing the unsteady glottal airflow are numerically solved in CFX using the finite volume method, while the finite element method (FEM) is used to solve the tissue deformation equations in the mechanical solver. Two data transfers are considered for the two-way FSI analysis: the first one is responsible to transfer the fluid force data from the CFD solver at the common airflow–VF tissue interface region to the fluid–solid interface region of the Mechanical solver. The second data transfer is for transferring the incremental displacement at fluid solid interface region in the mechanical system to the air domain as mesh displacement in the CFD system. The ANSYS System Coupling component controls the execution and convergence of fluid and solid simulations and solves fluid flow field and solid structure domains separately and in a sequential fashion. In particular, the module first performs the airflow simulations in ANSYS CFX to the target residual convergence to obtain the nonuniform pressure distribution and flow velocity in the larynx and provide the fluid forces exerted on the FSI interface. The obtained pressure load is then transferred to ANSYS mechanical to solve for corresponding deformations of VF tissue and obtain a converged solution of the displacement data. The obtained displacements are transferred back to the CFX module, and the problem is solved for convergence. The process is iterated for the number of defined coupling iterations till the completion of all the coupling steps. The maximum number of iterations per time step was set to 10 to ensure convergence within each time step. The coupling iterations are repeated until the convergence is reached, i.e., the interactions between fluid and structural components are converged. In System Coupling, a maximum root mean square residual of 0.01 must be reached to ensure the convergence of the solution of the coupled algorithm. These stagger iterations are reinitialized with the new deformed mesh that occurs from the solid deformation. After achieving the convergence of System Coupling, the glottal fluid flow forces of the FSI boundary are imported into the second CFX component, which independently solves the interstitial flow equations at each time step.

For the glottal flow and pressure analysis, the CFX module performs mesh deformation by solving a displacement diffusion equation to determine the mesh displacements throughout the remaining volume of the mesh. In particular, a locally increasing mesh stiffness (that is, the diffusivity for the mesh displacement equation) approach is used to provide additional control over the resulting mesh distribution and to avoid mesh folding and maintain an acceptable mesh quality in regions of large deformation. The existing solution will be interpolated to the updated mesh using tools provided with the ANSYS-CFX product. The prescribed coupling conditions for the continuity of velocities and the continuity of fluid and VF tissue stresses in the normal direction are controlled by the System Coupling data transfers. Data exchange is performed at the start of every coupling iteration within a coupling step. This process is handled with the Conservative Profile Preserving data transfer algorithm, and General Grid Interface (GGI) mapping algorithm is used to generate mapping weights when transferring forces. Moreover, the Profile Preserving algorithm and the Bucket Surface mapping algorithm are used to generate mapping weights when transferring displacement. The coupling box also manages incompatible meshes at the fluid–structure interface. All the simulations for the porous flow and pore pressure analysis are obtained using a fixed mesh algorithm, which somewhat eliminates the effect of the elastic property of the structure in filtration flow calculation. However, it should be noted that the fixed domain simplification has only been used for the porous flow analysis, and not for the FSI between the glottal flow and VF tissue.

The convergence criteria of the ANSYS CFX are chosen to 10⁻⁵ for the normalized residuals of the global linear system of equations for the mass and momentum. The convective term in the Navier–Stokes equations is linearized by the Picard iterations (equivalent to a fictitious time stepping method) namely, $v n + 1. ∇ v n + 1≈ v n. ∇ v n + 1$⁠, and the pressure variable is evaluated at the same nodes of the velocity. The Laplace operator in the momentum equations of the airflow domain is approximated by a centered scheme. A second-order upwind discretization scheme is applied for spatial discretization of the convective divergence terms, while an implicit second-order discretization scheme is used for the temporal approximation. The system is then solved using an algebraic multigrid method exploiting incomplete LU factorization as smoother. The DNS method for laminar flow is employed. The Reynolds number for the baseline case with $P i n$ = 1 kPa is estimated to be 1900 by considering the peak velocity of 38 m/s at the glottis and the maximum glottal gap of $d G max≅$ 0.8 mm during the VF opening phase of each vibrational cycle as the characteristic length.⁸⁴ This value is within the range for laryngeal airflow Reynolds numbers in previous computational⁸⁵ and experimental⁸⁶ studies.

The larynx model is discretized into unstructured linear tetrahedral elements within the airflow and VF regions in order to perform the simulations. The computational grid is selected after a mesh refinement study. This is accomplished by increasing the number of elements and conducting the simulations with refined grid density until the re-calculated model predictions agreed well with that from the original mesh. The resolution is considered sufficient for the VF region when the maximum VF displacement does not change appreciably (relative error less than 5%), and displacement analysis with coarser and finer meshes show negligible error when monitoring and comparing the difference between maximum values of glottal gap width and the total displacement patterns at some examined points. The fluid flow solution is also considered mesh-independent for an error lower than 5% in terms of velocity and pressure in the glottis. The final mesh is split into 1.2 × 10⁶ tetrahedral elements within the larynx and VF regions. Simulations are performed on a non-uniform grid which provides a higher resolution in the intraglottal region where high flow gradients are expected, as well as finer meshes in the groove areas and near the FSI interface of the VFs domain. A coarser mesh is used at the subglottal and supraglottal regions, decreasing the computation time and required resources. In most regions, the element size is 0.3–0.5 mm. Moreover, different size time steps are adopted to comply with the Cauchy–Freidrich–Lewey (CFL) numerical stability constraint and to provide adequate temporal resolution. A time step size of Δt = 10⁻⁵ s is applied on both the flow solver and the solid solver, which provides a good temporal resolution and an acceptable maximum CFL number with the optimal balance between computational efficiency and accuracy. Simulations are carried out whereby cases reached a well-defined state, such that the flow field is fully developed, and the VFs vibrate periodically.

Visualization of intraglottal aerodynamics by means of the velocity and pressure contours and the flow pattern using velocity streamlines in the investigated larynx model are provided in Fig. 3. The contours of the velocity and pressure, as well as the surface streamlines are obtained at the moment of glottal opening, in both the midcoronal and midsagittal planes.

The airflow is driven from the inlet to the outlet (left to right) by a pressure drop between the subglottal and supraglottal pressures and interacts with the permeable VFs to initiate the flow-induced vibrations. The spatial variation of air velocity in the larynx shows that the narrowing of the area for flow in the glottic region leads to a pressure drop within the glottal region [Fig. 3(b)] and therefore, an increase in intraglottal air velocity with the maximum value appeared between the folds [Fig. 3(a)]. More precisely, the velocity accelerates at the glottis where the diameter of the air lumen is smaller, while a drastic pressure gradient is observed where the velocity is the highest. The velocity streamlines indicate helical flow patterns and recirculation in the supraglottal region. The uniform streamlines at the inlet reach a maximum value between the folds at the glottis, forming a glottal jet, and then recirculate before exiting the outlet [Fig. 3(c)].

The results show that the model can reproduce the aerodynamic data through the larynx with a good approximant. In particular, the obtained pressure drops across the glottis and flow velocity fields are consistent with published modeling data^48,87,88 and experimental recording of glottal airflow.⁸⁰ 

Analysis of flow-induced vibration and the corresponding flow circulations in the VFs are explained in this section. Figure 4 features the deformation contours of the VF model at the AP mid surface, along with a view of the von Mises stress distribution and their corresponding time variations. Particularly, the profiles of total displacement pattern at the mid-coronal plane at five-time instants during one vibration cycle are displayed in Fig. 4(a). These spatiotemporal VF tissue deformations are obtained from opening to fully opened, then toward a closed position. Views of the total displacement demonstrate that the maximum deformations are at the medial VF surface. The subglottal pressure produces a pushing force in the flow direction on the VFs as well as a compressive force in the transverse direction, causing both the inferior and superior surfaces of the folds to abduct, and the glottis begins to open. The opening arrives at the maximum level, then adducting starts from the inferior area followed by the superior area and the VFs begin closing. Thus, the laryngeal airflow through the larynx causes the glottis opening and closing cycle to be repeated and the VFs to self-oscillate in a rhythmic vibration. This action of opening and closing corresponds to one phonation cycle, which causes the VF tissue to experience mechanical stresses.

During a phonation cycle [Fig. 4(a)], the $T 3$ time instant corresponds to the position of the VF when the glottis fully opens and deforms toward the superior direction (i.e., glottal gap has its maximum value). At this specific time point, the VF tissue von Mises stress (⁠ $σ V M$⁠) distribution contour on the midcoronal cross-section is shown in Fig. 4(b). The von Mises stress, also known as equivalent stress, is a measure of principal stresses denote by $S 1$⁠, $S 2$⁠, and $S 3$ and is calculated by $σ V M= 1 2 ( S 1 − S 2 ) 2 + S 3 − S 2 2 + ( S 1 − S 3 ) 2$⁠. It is observed that at the glottal opening the equivalent stress reaches the maximum value along the medial surface and a minimum value along the lateral surface.

Additionally, the time variations of the maximum total displacement and von Mises stress are plotted in Fig. 4(c). The phonation frequency is evaluated by plotting the variations of maximum VF displacement vs time [Fig. 4(c)] and indicates value of 120.1 Hz for the model, which agreed reasonably with the literature⁸⁹ for normal vocal range oscillation frequency. It can also be observed that for each phonation cycle, the VF vibration has reached a stationary state wherein the maximum displacement and therefore maximum glottal gap is nearly constant from cycle-to-cycle [Fig. 4(c)]. In particular, the vibration amplitude is approximately the same with less than 1 mm (0.63 mm) maximum deformation, highlighting the rhythmic pattern of phonation. The predicted vibration pattern mimics previous modeling⁶² and excised results.⁹⁰ Furthermore, the obtained displacement magnitudes are comparable with ones reported physiologically,⁹¹ and are also matched reasonably well with previous contours of VF displacement in models with similar lung pressure.^12,48,76 The variation of the equivalent stress mirrors that of the deformation with the maximum stress occurring at the peak deformation of each phonation cycle.

In addition to VF deformation and von Mises stress distribution contours, the interstitial fluid velocity and pore pressure are obtained. As the airflow causes the folds to self-oscillate, fluid circulates within the permeable VF tissue, and creates a filtration velocity and pressure. Figure 5 shows these interporous velocity vectors and pressure at the peak glottal flow, for the case with the porous properties listed in Table I. The velocity vectors obtained on the mid-coronal plane illustrate the circulation of liquid inside the tissue and the interporous liquid movement is represented as scaled vectors for ease of visualization [Fig. 5(a)]. Fluid circulation pattern in the VFs is also evaluated in three dimensions, distributed throughout the geometry. Color of the arrows represent the velocity magnitude, with red indicating the location of larger interstitial fluid flow. Rich spatial variation of the interstitial fluid velocity has been found [Fig. 5(a)]. Specifically, in the vertical direction, the fluid is transported from the inferior side to the superior side which conforms with inferior–superior redistribution observed in Ref. 35 during opening of VF. Results also indicate the fluid movement in the tissue is toward both the medial and lateral directions. The velocity reaches its highest magnitude across the VF medial surface where the fluid recirculates in the tissue, corresponding to the area where the von Mises stress is also the highest.

Furthermore, Fig. 5(b) shows the contour of pore pressure distribution depicted using a color gradient on both the coronal and transverse planes through the VF domain. The color bar represents the pressure magnitude with lower pressure values represented by colors toward the blue end of the spectrum, while red indicates higher values. Pore pressure is indicative of interporous fluid accumulation within the VF; a negative value indicates that the fluid is squeezed out of the region, whereas a positive pore pressure means that the fluid is accumulated. It is discovered that the inferior and medial surfaces experience the greatest pore pressure variations. The pore pressure gradient peaks at the AP midline because the fluid within the porous tissue moves toward the point of highest vibratory amplitude, consistent with previous observations in Refs. 44 and 92. The von Mises stress distribution follows a similar trend as the pore-pressure, increasing at the AP midline. The maximum pore pressure and filtration velocity within the VF agree well with previous studies in Refs. 43 and 44.

In this section, to investigate the importance of using a poroelastic model for the VF tissue, parametric studies focused on the impacts of lung pressure (⁠ $P i n$⁠) which provides the energy for VFs vibration, as well as the tissue permeability coefficient (⁠ $K V F$⁠) are performed. The objective is to quantify the resulting filtration flow and examine the range of model parameters that generate noticeable effects on the pore pressure and velocity. Four sets of lung pressure values and six sets of VF permeability coefficient are considered.

To investigate the effects of lung pressure and quantify the resulting liquid circulation, the subglottal inlet pressure varied for phonation loads of $P i n=1$⁠, 1.5, 2, and 2.5 kPa, while the VF permeability coefficient is kept constant at the value of $K V F= 10 − 13$ m², reported as the physiological permeability value (Table II). Using different lung pressures mentioned above, the FSI simulations between glottal airflow and VFs are performed. The chosen values of subglottal pressure cover the range of normal phonation conditions⁹³ and are comparable in magnitude to physiologically measured pressures used in previous models to achieve similar oscillation.^91,94 To determine the effect of fold permeability, the coefficient is varied in the range of $K V F$ = 10⁻¹⁵–10⁻⁷ m², with exponentially increasing values of $K V F= 10 − 15$⁠, $10 − 13$⁠, $10 − 11$⁠, $10 − 9$⁠, $10 − 8,$ and $10 − 7$ m², chosen according to VF permeability in previous studies,^35,44,78 while $P i n$ is kept constant at the previous specified value of 1 kPa. The coupled system is solved for each case to provide a prediction of the interstitial fluid flow within the tissue.

In Figs. 6–8, the vibration-induced flow circulation of VFs over a range of phonation conditions are analyzed. The filtration velocity streamlines, pore pressure, and tissue deformation on a surface in the AP midline for each different $P i n$ value at the fully opened instant are extracted and displayed in Fig. 6. In particular, Figs. 6(a) and 6(b) shows the comparison of the interstitial fluid velocity surface streamlines and pore pressure contours, while Fig. 6(c) illustrates the displacement distribution contours obtained in four different $P i n$ regimes to evaluate different phonation conditions. The maximum filtration velocity, pore pressure, and VF displacement are found at the medial surface in all the cases (Fig. 6). We also observe that both filtration velocity and pore pressure increase particularly along the medial surface of the VF with increasing lung pressure, while no change is observed across the lateral boundary. As more pressure is applied at the inlet of the laryngeal domain, the faster the fluid circulates within the VF tissue and a greater pressure gradient is obtained. Additionally, comparison of displacement distribution contours for all the cases indicates reduced levels of VF displacements as $P i n$ decreases, while the corresponding deformations of the VF in response to fluid induced vibrations increase significantly when the lung pressure is increased. The values of the total displacements of the VFs in millimeters are indicated by different scales of the color contours.

Furthermore, to observe the relationship between subglottal pressure, vibration amplitude, and subsequent liquid flow and stresses, the maximum VF displacement, maximum pore pressure, maximum von Mises stress, and maximum filtration velocity are computed and compared for each $P i n$ value. The relationships are displayed in Figs. 7 and 8. The maximum values during vibration are obtained from Fig. 6 when the folds are fully opened, and the variables $U max$⁠, maximum $σ V M$⁠, and maximum $P p$ for cases with varying inlet pressure are plotted (Fig. 7).

The maximum displacement at the AP midline of the fully opened VF at 1, 1.5, 2, and 2.5 kPa inlet pressures was at approximately 0.63, 0.97, 1.32, and 1.67 mm, respectively. Hence, as $P i n$ raised from 1 to 2.5 kPa, the maximum tissue displacement and therefore the vibration amplitude increases notably as expected and the vibration magnitude shows a linear relationship with inlet pressure (Fig. 7). A similar pattern is detected for the change in the maximum von Mises stress and pore pressure, revealing that pore pressure, tissue displacement, and von Mises stress all increase linearly with increasing the lung pressure, as shown in Fig. 7.

The variation of maximum and average filtration velocity in the VFs for cases with varying inlet lung pressures are illustrated in Fig. 8. The bars reflect the maximum filtration flow calculated at each subglottal pressure value and the line represents the change in the $q ¯ p$ trend. The left axis corresponds to the maximum velocity while the right axis depicts the average filtration velocity. Comparison of simulation results for maximum filtration velocity and volumetric average of filtration flow in the VF under varying $P i n$ pressure conditions show an increase in both. The volumetric average velocity shows a linear relationship with inlet pressure, while the maximum flow is detected when the lung pressure is highest.

We also compared the maximum filtration velocity for a range of tissue permeability values to quantify the extent that permeability coefficient affects the flow in the VF tissue. The results investigating the influence of $K V F$ are shown in Figs. 9 and 10. Figure 9 displays the filtration velocity streamlines in the porous VF tissue for the six cases with varying permeability values, and the legend is adjusted to show the variation with different ranges of velocity. As expected, the maximum filtration velocity during phonation increases with an increase in permeability coefficient and reaches the highest magnitude in the case where the $K V F$ is the highest which is shown in red (Fig. 9). Slightly different patterns for streamlines are noted for various permeability values. Comparing the results of pore pressure within the VF tissue (not visualized here) showed no significant changes in relation to VF permeability.

Additionally, the maximum filtration velocity as well as the volumetric average of filtration velocity over the VF domain calculated using Eq. (6) are compared for cases with different VF permeability coefficient (Fig. 10). The bars report the peak filtration flow, and the line represents the change in the average value, i.e., $q ¯ p$ trend. Results indicate that both the maximum filtration velocity and the volumetric average of filtration velocity vary exponentially with the change in permeability value. Namely, Fig. 10 reveals a positive nonlinear relationship between the two parameters: increasing the permeability increases the maximum and average of filtration velocity, the greatest peak filtration velocity value appears at the case where the permeability is the highest and decreases significantly as permeability decreases. Moreover, the cases with large permeability yield a significantly higher filtration flow within the entire region of the VF. Thus, the tissue permeability affects the filtration velocity distribution and its maximum value noticeably (Figs. 9 and 10) and plays a role in tissue perfusion during phonation.

The porous media plays an important role in oxygen supply and nutrient transport of soft tissue. Particularly, fluid circulation and the corresponding systemic hydration play a key chemo-mechanical role in the function of VFs, which is not fully understood. This fluid circulation and the subsequent rise in liquid pressure have been previously observed in a mechanical model of a vibrating capillary⁹² and is thought to contribute to benign VF injury development. Moreover, vibration of small blood vessels within the human VFs has been shown to vary perfusion flow rate,⁹⁵ which in turn influences processes, such as nutrient and waste transport.

The current study aims to account for the poroelasticity of the tissue and merge that with the interaction between glottal airflow and VF dynamics, in a FPSI computational framework which allows for modeling the interstitial fluid dynamics and perfusion in vibrating VFs. In all the previous FSI simulations of phonation, the fold is treated as an impermeable single-phase structure due to significant complexity of incorporating a biphasic tissue. In particular, self-oscillation results from coupling of airflow to VF movement and provides a physiologically relevant way to explore the effects of airflow aerodynamics and vary boundary loads not only with time but also across the surface; however, it is computationally costly. Because of the high computational demands of modeling the interaction of the airflow, the poroelastic solid skeleton, and the interporous liquid, boundary loads in all the previous porous VF models have been applied directly, which can sacrifice the qualitative outcome. The developed FPSI computational model in this manuscript provides valuable insight into the salient behavior of interstitial fluid circulation which is a critical metric of the tissue's functional state, and a prediction of filtration velocity within the porous VFs at different phonation conditions and porous features of the tissue. Particularly, the extent that VF permeability variation (⁠ $K V F$⁠) and phonation conditions (⁠ $P i n$⁠) affect the filtration velocity is quantified to highlight the importance of including poroelasticity in phonation models.

The simulation results provided qualitative and quantitative prediction of filtration flow patterns and pressure distribution in the VF tissue. Velocity streamlines as well as both pressure and velocity contours of the glottal airflow through the larynx are also obtained and analyzed. The interstitial fluid flow is found to be heterogeneous with its maximum value at the AP midline of the folds, as well as upward movement in the inferior–superior direction similar to previous observations by Ref. 44. It has been also noted that this filtration velocity is strongly dependent on both the phonation condition as reported previously in Ref. 96, as well as tissue permeability coefficient. A linear relationship is found between the vibration amplitude, pore pressure, and filtration flow with subglottal pressure, comparable with Ref. 44, which is likely due to the fact that the presented model is perfectly elastic. The obtained linear relationship between the vibration magnitude and the subglottal pressure is also consistent with previous numerical findings⁹⁷ for the similar range of lung pressure used in the current study. Although findings revealed that the maximum filtration velocity decreases with reducing the subglottal pressure, insignificant differences in flow pattern [indicated by the streamlines in Fig. 6(a)] are noted when comparing the cases with varying lung pressure.

Furthermore, VF porous properties significantly impact the perfusion and filtration velocity within the tissue and significant differences are noted when comparing cases with different permeability coefficient values. Namely, the maximum and average filtration velocities varied greatly among the cases with different permeabilities (Figs. 9 and 10). The relationship between filtration velocity and permeability is found to be nonlinear. Larger recirculation vortex zones are observed for cases with higher VF permeability that shifted the flow fields toward the medial surface of the VFs. However, no significant differences appeared for the total wall displacement and pore pressure when the results between the cases with varying $K V F$ are compared.

Like most studies attempting to simulate an in vivo condition, this investigation has limitations that need to be considered while interpreting the results. While similar to the approach taken in past VF studies an idealized model of the larynx is employed, the inclusion of patient-specific geometries constructed from medical imaging is necessary to establish models that predict the phonation conditions that are benign to voice health. Another limitation is the imposed artificial gap between the two folds, which allows for flow leakage during a glottal closure. Although the presence of structural contact and a complete glottal closure is likely to affect the stresses in VF, the purpose of this study is to investigate the effect of including a porous tissue on filtration velocity and flow distribution, in particular when VF is opened. Moreover, the fluid exchange condition through the VF lateral boundary is another parametric factor that influences the flow pattern within a VF, as it decides the amount of fluid entering and exiting the porous media. Also, the predicted filtration fluid circulation is entirely driven by the interactions of the airflow on VFs. Namely, the glottal flow forces on the tissue squeezes the fluid and causes its movement. To make this computational model functional in clinical setting, the prevalence of VF vibration in filtration fluid dynamics as well as the frequency analysis of tissue vibration and the resulting acoustic measures need to be incorporated within the framework to comprehensively investigate the relationship between VF porous properties and the vibration mode, which can be the subject of future research.

Additionally, VF permeability coefficient and porosity lack experimental measurements. The value of these poroelastic parameters could influence the flow property of fluid and smaller permeability values downgrade VF into a classic elastic model. Moreover, the simulations are performed for a limited number of values for each parameter to represent the variations of phonation conditions and porous features of tissue. Although the results support the important role of permeability on VF biomechanics, it would be ideal to consider a wider range of parameters to reach stronger conclusions. Furthermore, tissue permeability is assumed uniform throughout the VF domain due to the limitations in the experimental data, while it has been shown that the vast majority of fibers in VF tissue are aligned longitudinally, and the permeability of liquid parallel to the fibers is larger than in the transverse directions. In addition, the VF tissue skeleton is assumed to be homogeneous and isotropic, with a linear stress–strain relationship, while experimental evidence suggests that the VF is an anisotropic material. A linear elastic approximation is valid for the observed deformations and has worked well for previous models. However, this can affect the realistic dynamics of VFs during the closing phase. The degree of error introduced by these assumptions will be investigated in future studies. However, it should be noted that despite the limitations that these assumptions could have imposed, results in this study are consistent with previous findings as explained in Sec. III. Incorporation of aforementioned enhancements would not alter the current conclusion that VF permeability and subglottal pressure affect the filtration flow in the tissue and the values of comparative results are still valid, considering that the computational cost of the model should be balanced with accuracy.

Another limitation is assuming the laminar flow for laryngeal models. In this study, airflow Reynolds number in larynx for various cases of inlet pressure (for the simulation at the low-pressure and high-pressure conditions) is observed to be in the range of $≈$ 1900–4000 for the glottal gap and maximum velocity measured in each case, consistent with the obtained range of 1000–5000 in previous research.⁸⁶ However, comprehensive investigations have revealed complex flow patterns and morphologies near the glottal region⁹⁸ characterized by features such as transition to turbulence and the corresponding vortex structures.⁹⁹ The numerical computation of this highly unsteady airflow with massive separation is the focus of future studies. Furthermore, although encouraging results are obtained, modeling the VFs as a three-layered transversely isotropic porous tissue, representing the body, ligament, and cover will be considered as a future direction to improve the modeling of material properties. Future work also extends this model to predict oxygen distribution in the porous medium of VF tissue considering diffusion and convection, similar to previous coupled mass transport–fluid–porous structure simulations,^38,100,101 and explore the flow within the folds and its relation to hypoxia,^102,103 which can also provide additional information on VF hydration.^39,104

Although this study had certain computational simplifications, it presents the first biphasic FSI model that uses a realistic VF geometry and characterizes perfusion through a vibrating tissue. Our observations suggest that flow to the VF increases during phonation, which is in agreement with¹⁰⁵ that directly measured the blood flow to the VFs using microsphere surface technique. Moreover, while reducing the lung pressure resulted in the weakening of perfusion flow in the tissue, insignificant differences in flow recirculation patterns are noted by varying subglottal pressure. The poroelastic theory provides information about interstitial fluid dynamics in the vibrating VFs which cannot be predicted by the often-used elastic theory. Further investigation of these mechanical factors may lend insight into the mechanism of benign lesion formation,⁹² optimizing the design of laryngeal prostheses¹⁰⁶ and disease treatment, and exploring hemodynamics in the tissue as well as the flow combined with drug delivery.¹⁰⁷ 

This study was supported by CBET 2138225 from the National Science Foundation.

The authors have no conflicts to disclose.

Isabella McCollum: Data curation (equal); Formal analysis (lead); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Durwash Badr: Data curation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Alexis Throop: Data curation (equal); Methodology (supporting); Visualization (equal). Rana Zakerzadeh: Conceptualization (lead); Funding acquisition (lead); Methodology (lead); Project administration (lead); Supervision (lead); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.