An immersed-boundary method for flow-structure interaction in biological systems with application to phonation

An Integrated Experimental-Computational Study of Vocal Fold Vibration in Type I Thyroplasty

Subject-specific computational modeling of vocal fold (VF) vibration was integrated with an ex vivo animal experiment of type 1 thyroplasty to study the effect of the implant on the vocal fold vibration. In the experiment, a rabbit larynx was used to simulate type 1 thyroplasty, where one side of the vocal fold was medialized with a trans-muscular suture while the other side was medialized with a silastic implant. Vocal fold vibration was then achieved by flowing air through the larynx and was filmed with a high-speed camera. The three-dimensional computational model was built upon the pre-operative scan of the laryngeal anatomy. This subject-specific model was used to simulate the vocal fold medialization and then the fluid-structure interaction (FSI) of the vocal fold. Model validation was done by comparing the vocal fold displacement with postoperative scan (for medialization), and by comparing the vibratory characteristics with the high-speed images (for vibration). These comparisons showed the computational model successfully captured the effect of the implant and thus has the potential for presurgical planning.

The effect of swelling on vocal fold kinematics and dynamics

Swelling in the vocal folds is caused by the local accumulation of fluid, and has been implicated as a phase in the development of phonotraumatic vocal hyperfunction and related structural pathologies, such as vocal fold nodules. It has been posited that small degrees of swelling may be protective, but large amounts may lead to a vicious cycle wherein the engorged folds lead to conditions that promote further swelling, leading to pathologies. As a first effort to explore the mechanics of vocal fold swelling and its potential role in the etiology of voice disorders, this study employs a finite-element model with swelling confined to the superficial lamina propria, which changes the volume, mass, and stiffness of the cover layer. The impacts of swelling on a number of vocal fold kinematic and damage measures, including von Mises stress, internal viscous dissipation, and collision pressure, are presented. Swelling has small but consistent effects on voice outputs, including a reduction in fundamental frequency with increasing swelling (10 Hz at 30 % swelling). Average von Mises stress decreases slightly for small degrees of swelling but increases at large magnitudes, consistent with expectations for a vicious cycle. Both viscous dissipation and collision pressure consistently increase with the magnitude of swelling. This first effort at modeling the impact of swelling on vocal fold kinematics, kinetics, and damage measures highlights the complexity with which phonotrauma can influence performance metrics. Further identification and exploration of salient candidate measures of damage and refined studies coupling swelling with local phonotrauma are expected to shed further light on the etiological pathways of phonotraumatic vocal hyperfunction.

The Influence of Fiber Orientation of the Conus Elasticus in Vocal Fold Modeling

While the conus elasticus is generally considered a part of continuation of the vocal ligament, histological studies have revealed different fiber orientations that fibers are primarily aligned in the superior-inferior direction in the conus elasticus and in the anterior-posterior direction in the vocal ligament. In this work, two continuum vocal fold models are constructed with two different fiber orientations in the conus elasticus: the superior-inferior direction and the anterior-posterior direction. Flow-structure interaction simulations are conducted at different subglottal pressures to investigate the effects of fiber orientation in the conus elasticus on vocal fold vibrations, aerodynamic and acoustic measures of voice production. The results reveal that including the realistic fiber orientation (superior-inferior) in the conus elasticus yields smaller stiffness and larger deflection in the coronal plane at the junction of the conus elasticus and ligament and subsequently leads to a greater vibration amplitude and larger mucosal wave amplitude of the vocal fold. The smaller coronal-plane stiffness also causes a larger peak flow rate and higher skewing quotient. Furthermore, the voice generated by the vocal fold model with a realistic conus elasticus has a lower fundamental frequency, smaller first harmonic amplitude, and smaller spectral slope.

Towards modeling of phonation and its recovery in unilateral vocal fold paralysis by fluid-structure interaction

Introduction

Vocal folds are responsible for sound generation. In unilateral vocal fold paralysis (UVFP), the recurrent laryngeal nerve, which controls the vocal folds, is paralyzed. Medialization laryngoplasty is a surgery in which an implant is inserted to push the paralyzed vocal fold to the centerline to recover phonation.

Methods

Here, a numerical simulation is used to calculate flow-related parameters to give insight into what happens in healthy and treated(implanted) vocal folds and their enhancement. In the present work, airflow over vocal folds is modeled considering fluid-structure interaction (FSI) and varying inlet pressure. The governing equations are discretized for fluid and solid domains and solved using the Galerkin finite element method. The boundary conditions for healthy and unilaterally paralyzed vocal folds were imposed to agree with real cases behavior.

Results

The results showed the effectiveness of medialization laryngoplasty in treating unilateral vocal fold paralysis concerning healthy vocal folds.

Conclusion

This simulation provided a better insight into treatment results for patient-specific cases.

A Deep Learning-Based Generalized Empirical Flow Model of Glottal Flow During Normal Phonation

This paper proposes a deep learning-based generalized empirical flow model (EFM) that can provide a fast and accurate prediction of the glottal flow during normal phonation. The approach is based on the assumption that the vibration of the vocal folds can be represented by a universal kinematics equation (UKE), which is used to generate a glottal shape library. For each shape in the library, the ground truth values of the flow rate and pressure distribution are obtained from the high-fidelity Navier-Stokes (N-S) solution. A fully connected deep neural network (DNN) is then trained to build the empirical mapping between the shapes and the flow rate and pressure distributions. The obtained DNN-based EFM is coupled with a finite element method (FEM)-based solid dynamics solver for fluid-structure-interaction (FSI) simulation of phonation. The EFM is evaluated by comparing the N-S solutions in both static glottal shapes and FSI simulations. The results demonstrate a good prediction performance in accuracy and efficiency.

An Improved Glottal Flow Model Based on Seq2Seq LSTM for Simulation of Vocal Fold Vibration

OBJECTIVES

An improved data-driven glottal flow model for fluid-structure interaction (FSI) simulation of the vocal fold vibration is proposed in this paper. This model aims to improve the prediction performance of the previously developed deep neural network (DNN) based empirical flow model (EFM)1 on accuracy and efficiency.

METHODS

A Seq2Seq long short-term memory (LSTM) network is employed in the present model to infer the flow rate and pressure distribution from the subglottal pressure and cross-section area distribution of the glottis. The training data is collected from the generalized glottal shape library generated in Zhang et al.1 RESULTS AND CONCLUSIONS: Compared to the EFM, the present model not only discards the time-consuming optimization process, but also drastically reduces the errors, therefore the prediction performance can be greatly improved. The present model is evaluated by coupling with a solid dynamics solver for FSI simulation, and the results demonstrate a great improvement on accuracy and efficiency.

Subject-Specific Computational Fluid-Structure Interaction Modeling of Rabbit Vocal Fold Vibration

A full three-dimensional (3D) fluid-structure interaction (FSI) study of subject-specific vocal fold vibration is carried out based on the previously reconstructed vocal fold models of rabbit larynges. Our primary focuses are the vibration characteristics of the vocal fold, the unsteady 3D flow field, and comparison with a recently developed 1D glottal flow model that incorporates machine learning. The 3D FSI model applies strong coupling between the finite-element model for the vocal fold tissue and the incompressible Navier-Stokes equation for the flow. Five different samples of the rabbit larynx, reconstructed from the magnetic resonance imaging (MRI) scans after the in vivo phonation experiments, are used in the FSI simulation. These samples have distinct geometries and a different inlet pressure measured in the experiment. Furthermore, the material properties of the vocal fold tissue were determined previously for each individual sample. The results demonstrate that the vibration and the intraglottal pressure from the 3D flow simulation agree well with those from the 1D flow model based simulation. Further 3D analyses show that the inferior and supraglottal geometries play significant roles in the FSI process. Similarity of the flow pattern with the human vocal fold is discussed. This study supports the effective usage of rabbit larynges to understand human phonation and will help guide our future computational studies that address vocal fold disorders.

Computational Modeling of Voice Production Using Excised Canine Larynx

A combined experimental-numerical work was conducted to comprehensively validate a subject-specific continuum model of voice production in larynx using excised canine laryngeal experiments. The computational model is a coupling of the Navier-Stokes equations for glottal flow dynamics and a finite element model of vocal fold dynamics. The numerical simulations employed a cover-body vocal fold structure with the geometry reconstructed from magnetic resonance imaging scans and the material properties determined through an optimization-based inverse process of experimental indentation measurement. The results showed that the simulations predicted key features of the dynamics observed in the experiments, including the skewing of the glottal flow waveform, mucosal wave propagation, continuous increase of the divergent angle and intraglottal swirl strength during glottal closing, and flow recirculation between glottal jet and vocal fold. The simulations also predicted the increase of the divergent angle, glottal jet speed, and intraglottal flow swirl strength with the subglottal pressure, same as in the experiments. Quantitatively, the simulations over-predicted the frequency and jet speed and under-predicted the flow rate and divergent angle for the larynx under study. The limitations of the model and their implications were discussed.

Automatic triangulated mesh generation of pulmonary airways from segmented lung 3DCTs for computational fluid dynamics

PURPOSE

Computational fluid dynamics (CFD) of lung airflow during normal and pathophysiological breathing provides insight into regional pulmonary ventilation. By integrating CFD methods with 4D lung imaging workflows, regions of normal pulmonary function can be spared during treatment planning. To facilitate the use of CFD simulations in a clinical setup, a robust, automated, and CFD-compliant airway mesh generation technique is necessary.

METHODS

We define a CFD-compliant airway mesh to be devoid of blockages of airflow and leaks in the airway path, both of which are caused by airway meshing errors that occur when using conventional meshing techniques. We present an algorithm to create a CFD-compliant airway mesh in an automated manner. Beginning with a medial skeleton of the airway segmentation, the branches were tracked, and 3D points at which bifurcations occur were identified. Airway branches and bifurcation features were isolated to allow for automated and careful meshing that considered their anatomical nature.

RESULTS

We present the meshing results from three state-of-the-art tools and compare them with the meshes generated by our algorithm. The results show that fully CFD-compliant meshes were automatically generated for an ideal geometry and patient-specific CT scans. Using an open-source smoothed-particle hydrodynamics CFD implementation, we compared the airflow using our approach and conventionally generated airway meshes.

CONCLUSION

Our meshing algorithm was able to successfully generate a CFD-compliant mesh from pre-segmented lung CT scans, providing an automatic meshing approach that enables interventional CFD simulations to guide lung procedures such as radiotherapy or lung volume reduction surgery.

Effects of Sulcus Vocalis Depth on Phonation in Three-Dimensional Fluid-Structure Interaction Laryngeal Models

Sulcus vocalis is an indentation parallel to the edge of vocal fold, which may extend into the cover and ligament layer of the vocal fold or deeper. The effects of sulcus vocalis depth d on phonation and the vocal cord vibrations are investigated in this study. The three-dimensional laryngeal models were established for healthy vocal folds (0 mm) and different types of sulcus vocalis with the typical depth of 1 mm, 2 mm, and 3 mm. These models with fluid-structure interaction (FSI) are computed numerically by sequential coupling method, which includes an immersed boundary method (IBM) for modelling the glottal airflow, a finite-element method (FEM) for modelling vocal fold tissue. The results show that a deeper sulcus vocalis in the cover layer decreases the vibrating frequency of vocal folds and expands the prephonatory glottal half-width which increases the phonation threshold pressure. The larger sulcus vocalis depth makes vocal folds difficult to vibrate and phonate. The effects of sulcus vocalis depth suggest that the feature such as phonation threshold pressure could assist in the detection of healthy vocal folds and different types of sulcus vocalis.

3D-FV-FE Aeroacoustic Larynx Model for Investigation of Functional Based Voice Disorders

For the clinical analysis of underlying mechanisms of voice disorders, we developed a numerical aeroacoustic larynx model, called simVoice, that mimics commonly observed functional laryngeal disorders as glottal insufficiency and vibrational left-right asymmetries. The model is a combination of the Finite Volume (FV) CFD solver Star-CCM+ and the Finite Element (FE) aeroacoustic solver CFS++. simVoice models turbulence using Large Eddy Simulations (LES) and the acoustic wave propagation with the perturbed convective wave equation (PCWE). Its geometry corresponds to a simplified larynx and a vocal tract model representing the vowel /a/. The oscillations of the vocal folds are externally driven. In total, 10 configurations with different degrees of functional-based disorders were simulated and analyzed. The energy transfer between the glottal airflow and the vocal folds decreases with an increasing glottal insufficiency and potentially reflects the higher effort during speech for patients being concerned. This loss of energy transfer may also have an essential influence on the quality of the sound signal as expressed by decreasing sound pressure level (SPL), Cepstral Peak Prominence (CPP), and Vocal Efficiency (VE). Asymmetry in the vocal fold oscillations also reduces the quality of the sound signal. However, simVoice confirmed previous clinical and experimental observations that a high level of glottal insufficiency worsens the acoustic signal quality more than oscillatory left-right asymmetry. Both symptoms in combination will further reduce the quality of the sound signal. In summary, simVoice allows for detailed analysis of the origins of disordered voice production and hence fosters the further understanding of laryngeal physiology, including occurring dependencies. A current walltime of 10 h/cycle is, with a prospective increase in computing power, auspicious for a future clinical use of simVoice.

A one-dimensional flow model enhanced by machine learning for simulation of vocal fold vibration

A one-dimensional (1D) unsteady and viscous flow model that is derived from the momentum and mass conservation equations is described, and to enhance this physics-based model, a machine learning approach is used to determine the unknown modeling parameters. Specifically, an idealized larynx model is constructed and ten cases of three-dimensional (3D) fluid-structure interaction (FSI) simulations are performed. The flow data are then extracted to train the 1D flow model using a sparse identification approach for nonlinear dynamical systems. As a result of training, we obtain the analytical expressions for the entrance effect and pressure loss in the glottis, which are then incorporated in the flow model to conveniently handle different glottal shapes due to vocal fold vibration. We apply the enhanced 1D flow model in the FSI simulation of both idealized vocal fold geometries and subject-specific anatomical geometries reconstructed from the magnetic resonance imaging images of rabbits' larynges. The 1D flow model is evaluated in both of these setups and shown to have robust performance. Therefore, it provides a fast simulation tool that is superior to the previous 1D models.

Immersed boundary method for multiphase transport phenomena

Multiphase flows with momentum, heat, and mass transfer exist widely in a variety of industrial applications. With the rapid development of numerical algorithms and computer capacity, advanced numerical simulation has become a promising tool in investigating multiphase transport problems. Immersed boundary (IB) method has recently emerged as such a popular interface capturing method for efficient simulations of multiphase flows, and significant achievements have been obtained. In this review, we attempt to give an overview of recent progresses on IB method for multiphase transport phenomena. Firstly, the governing equations, the basic ideas, and different boundary conditions for the IB methods are introduced. This is followed by numerical strategies, from which the IB methods are classified into two types, namely the artificial boundary method and the authentic boundary method. Discussions on the implementation of various boundary conditions at the interphase surface with momentum, heat, and mass transfer for different IB methods are then presented, together with a summary. Then, the state-of-the-art applications of IB methods to multiphase flows, including the isothermal flows, the heat transfer flows, and the mass transfer problems are outlined, with particular emphasis on the latter two topics. Finally, the conclusions and future challenges are identified.

Numerical study of droplet impact on a flexible substrate

Droplets interacting with deformable moving boundaries is ubiquitous. The flexible boundaries may dramatically affect the hydrodynamic behavior of droplets. A numerical method for simulating droplet impact on flexible substrates is developed. The effect of flexibility is investigated. To reduce the contact time and increase the remaining upward momentum in the flexible cases, the Weber number should be larger than a critical value. Moreover, the ratio of the natural frequency of the plate to that of the droplet F_{r} should approximately equal to the reciprocal of the contact time of droplets impact on the rigid surfaces (t_{ctr}) at the same We, e.g., F_{r}≈1/t_{ctr}. Only under this circumstance would the kinetic energy convert into the surface energy of the droplet and the elastic energy of the plate simultaneously, and vice versa. Moreover, based on a double spring model, we proposed scaling laws for the maximal deflection of the plate and spreading diameter of the drop. Finally, the droplet impact under different wettability is qualitatively studied. We found that the flexibility may contribute to the droplet bouncing at a smaller contact angle.

Quantification of the Intraglottal Pressure Induced by Flow Separation Vortices Using Large Eddy Simulation

The greatest rate of change in the glottal flow rate during phonation is a rapid decrease that occurs during the latter part of the glottal closing. Previous works showed that intraglottal flow separation vortices form in a divergent glottis, produce negative gauge pressures (below atmospheric) during closing. It is hypothesized here that flow separation vortices contribute to the rapid closing mechanism of the true vocal folds during phonation. Four idealized static models (M5) of the human larynx were investigated using large eddy simulation: 2 models featured parallel folds that did not enable flow separation in the glottis and 2 models involved a divergent glottis. The influence of the ventricular gap (narrow/wide) is evaluated. An unsteady pressure inlet representing a voicing cycle was applied to the sub-glottal region to mimic the time-varying glottal flow. Intraglottal vortex structures formed downstream of the separation point in a divergent glottis. Their existence caused a higher closing force that was applied onto the vocal folds. A narrow ventricular gap strengthens this effect. Strength of the intraglottal vortices increased with the maximum flow declination rate. Therefore, a more divergent shape of the glottis during glottal closing will be one of the main contributors to the voice quality.

A reduced-order flow model for vocal fold vibration: from idealized to subject-specific models

We present a reduced-order model for fluid-structure interaction (FSI) simulation of vocal fold vibration during phonation. This model couples the three-dimensional (3D) tissue mechanics and a one-dimensional (1D) flow model that is derived from the momentum and mass conservation equations for the glottal airflow. The effects of glottal entrance and pressure loss in the glottis are incorporated in the flow model. We consider both idealized vocal fold geometries and subject-specific anatomical geometries segmented from the MRI images of rabbits. For the idealized vocal fold geometries, we compare the simulation results from the 1D/3D hybrid FSI model with those from the full 3D FSI simulation based on an immersed-boundary method. For the subject-specific geometries, we incorporate previously estimated tissue properties for individual samples and compare the results with those from the high-speed imaging experiment of in vivo phonation. In both setups, the comparison shows good agreement in the vibration frequency, amplitude, phase delay, and deformation pattern of the vocal fold, which suggests potential application of the present approach for future patient-specific modeling.

Volume velocity in a canine larynx model using time‑resolved tomographic particle image velocimetry

In the classic source-filter theory, the source of sound is flow modulation. "Flow" is the flow rate (Q) and flow modulation is dQ/dt. Other investigators have argued, using theoretical, computational, and mechanical models of the larynx, that there are additional sources of sound. To determine the acoustic role of dQ/dt in a tissue model, Q needs to be accurately measured within a few millimeters of the glottal exit; however, no direct measures of Q currently exist. The goal of this study is to obtain this waveform in an excised canine larynx model using time-resolved tomographic particle image velocimetry. The flow rate data are captured simultaneously with acoustic measurements to determine relations with vocal characteristics. The results show that glottal waveform characteristics such as maximum flow declination rate are proportional to the subglottal pressure, fundamental frequency, and acoustic intensity. These findings are important as they use direct measurements of the volume flow at the glottal exit to validate some of the assumptions used in the source-filter theory. In addition, future work will address the accuracy of indirect clinical measurement techniques, such as the Rothenberg mask.

A Reduced-Order Flow Model for Fluid-Structure Interaction Simulation of Vocal Fold Vibration

We present a novel reduced-order glottal airflow model that can be coupled with the three-dimensional (3D) solid mechanics model of the vocal fold tissue to simulate the fluid-structure interaction (FSI) during voice production. This type of hybrid FSI models have potential applications in the estimation of the tissue properties that are unknown due to patient variations and/or neuromuscular activities. In this work, the flow is simplified to a one-dimensional (1D) momentum equation-based model incorporating the entrance effect and energy loss in the glottis. The performance of the flow model is assessed using a simplified yet 3D vocal fold configuration. We use the immersed-boundary method-based 3D FSI simulation as a benchmark to evaluate the momentum-based model as well as the Bernoulli-based 1D flow models. The results show that the new model has significantly better performance than the Bernoulli models in terms of prediction about the vocal fold vibration frequency, amplitude, and phase delay. Furthermore, the comparison results are consistent for different medial thicknesses of the vocal fold, subglottal pressures, and tissue material behaviors, indicating that the new model has better robustness than previous reduced-order models.

A Deep Neural Network Based Glottal Flow Model for Predicting Fluid-Structure Interactions during Voice Production

This paper proposes a machine-learning based reduced-order model that can provide fast and accurate prediction of the glottal flow during voice production. The model is based on the Bernoulli equation with a viscous loss term predicted by a deep neural network (DNN) model. The training data of the DNN model is a Navier-Stokes (N-S) equation-based three-dimensional simulation of glottal flows in various glottal shapes generated by a synthetic shape function, which can be obtained by superimposing the instantaneous modal displacements during vibration on the prephonatory geometry of the glottal shape. The input parameters of the DNN model are the geometric and flow parameters extracted from discretized cross sections of the glottal shapes and the output target is the corresponding flow resistance coefficient. With this trained DNN-Bernoulli model, the flow resistance coefficient as well as the flow rate and pressure distribution in any given glottal shape generated by the synthetic shape function can be predicted. The model is further coupled with a finite-element method based solid dynamics solver for simulating fluid-structure interactions (FSI). The prediction performance of the model for both static shape and FSI simulations is evaluated by comparing the solutions to those obtained by the Bernoulli and N-S model. The model shows a good prediction performance in accuracy and efficiency, suggesting a promise for future clinical use.

A Deep Neural Network Based Glottal Flow Model for Predicting Fluid-Structure Interactions during Voice Production

This paper proposes a machine-learning based reduced-order model that can provide fast and accurate prediction of the glottal flow during voice production. The model is based on the Bernoulli equation with a viscous loss term predicted by a deep neural network (DNN) model. The training data of the DNN model is a Navier-Stokes (N-S) equation-based three-dimensional simulation of glottal flows in various glottal shapes generated by a synthetic shape function, which can be obtained by superimposing the instantaneous modal displacements during vibration on the prephonatory geometry of the glottal shape. The input parameters of the DNN model are the geometric and flow parameters extracted from discretized cross sections of the glottal shapes and the output target is the corresponding flow resistance coefficient. With this trained DNN-Bernoulli model, the flow resistance coefficient as well as the flow rate and pressure distribution in any given glottal shape generated by the synthetic shape function can be predicted. The model is further coupled with a finite-element method based solid dynamics solver for simulating fluid-structure interactions (FSI). The prediction performance of the model for both static shape and FSI simulations is evaluated by comparing the solutions to those obtained by the Bernoulli and N-S model. The model shows a good prediction performance in accuracy and efficiency, suggesting a promise for future clinical use.

Three-Dimensional Numerical Analysis of Periciliary Liquid Layer: Ciliary Abnormalities in Respiratory Diseases

Human pulmonary epithelial cells are protected by two layers of fluid—the outer watery periciliary liquid layer (PCL) and the uppermost non-Newtonian mucus layer (ML). Aerosols and inhaled toxic particles are trapped by the ML which must then be removed swiftly to avoid adverse health implications. Epithelial cells are covered with cilia that beat rapidly within the PCL. Such ciliary motion drives the mucus transport. Although cilia can penetrate slightly inside the mucus to assist mucus movement, the motion of the underlying PCL layer within the airway surface liquid (ASL) is significant in mucus and pathogens transport. As such, a detailed parametric study of the influence of different abnormal cilia characteristics, such as low beating frequency, short length, abnormal beating pattern, reduced ciliary density, and epithelium patchiness due to missing cilia on the PCL transport, is carried out numerically. Such abnormalities are found in various chronic respiratory diseases. In addition, the shear stress at the epithelium is assessed due to the importance of shear stress on the epithelial function. Using the immersed boundary (IB) method combined with the finite-difference projection method, we found that the PCL, under standard healthy conditions, has net forward motion but that different diseased conditions decrease the forward motion of the PCL, as is expected based on clinical understanding.

Immersed Methods for Fluid-Structure Interaction

Fluid-structure interaction is ubiquitous in nature and occurs at all biological scales. Immersed methods provide mathematical and computational frameworks for modeling fluid-structure systems. These methods, which typically use an Eulerian description of the fluid and a Lagrangian description of the structure, can treat thin immersed boundaries and volumetric bodies, and they can model structures that are flexible or rigid or that move with prescribed deformational kinematics. Immersed formulations do not require body-fitted discretizations and thereby avoid the frequent grid regeneration that can otherwise be required for models involving large deformations and displacements. This article reviews immersed methods for both elastic structures and structures with prescribed kinematics. It considers formulations using integral operators to connect the Eulerian and Lagrangian frames and methods that directly apply jump conditions along fluid-structure interfaces. Benchmark problems demonstrate the effectiveness of these methods, and selected applications at Reynolds numbers up to approximately 20,000 highlight their impact in biological and biomedical modeling and simulation.

Finite Element Framework for Computational Fluid Dynamics in FEBio

The mechanics of biological fluids is an important topic in biomechanics, often requiring the use of computational tools to analyze problems with realistic geometries and material properties. This study describes the formulation and implementation of a finite element framework for computational fluid dynamics (CFD) in FEBio, a free software designed to meet the computational needs of the biomechanics and biophysics communities. This formulation models nearly incompressible flow with a compressible isothermal formulation that uses a physically realistic value for the fluid bulk modulus. It employs fluid velocity and dilatation as essential variables: The virtual work integral enforces the balance of linear momentum and the kinematic constraint between fluid velocity and dilatation, while fluid density varies with dilatation as prescribed by the axiom of mass balance. Using this approach, equal-order interpolations may be used for both essential variables over each element, contrary to traditional mixed formulations that must explicitly satisfy the inf-sup condition. The formulation accommodates Newtonian and non-Newtonian viscous responses as well as inviscid fluids. The efficiency of numerical solutions is enhanced using Broyden's quasi-Newton method. The results of finite element simulations were verified using well-documented benchmark problems as well as comparisons with other free and commercial codes. These analyses demonstrated that the novel formulation introduced in FEBio could successfully reproduce the results of other codes. The analogy between this CFD formulation and standard finite element formulations for solid mechanics makes it suitable for future extension to fluid-structure interactions (FSIs).

A Formulation for Fluid Structure-Interactions in FEBio Using Mixture Theory

Many physiological systems involve strong interactions between fluids and solids, posing a signicant challenge when modeling biomechanics. The objective of this study was to implement a fluid-structure interaction (FSI) solver in the free, open-source finite element code FEBio (febio.org), that combined the existing solid mechanics and rigid body dynamics solver with a recently-developed computational fluid dynamics (CFD) solver. A novel Galerkin-based finite element FSI formulation was introduced based on mixture theory, where the FSI domain was described as a mixture of fluid and solid constituents that have distinct motions. The mesh was defined on the solid domain, specialized to have zero mass, negligible stiffness and zero frictional interactions with the fluid, whereas the fluid was modeled as isothermal and compressible. The mixture framework provided the foundation for evaluating material time derivatives in a material frame for the solid and in a spatial frame for the fluid. Similar to our recently reported CFD solver, our FSI formulation did not require stabilization methods to achieve good convergence, producing a compact set of equations and code implementation. The code was successfully verified against benchmark problems and an analytical solution for squeeze-film lubrication. It was validated against experimental measurements of the flow rate in a peristaltic pump, and illustrated using non-Newtonian blood flow through a bifurcated carotid artery with a thick arterial wall. The successful formulation and implementation of this FSI solver enhances the multiphysics modeling capabilities in FEBio relevant to the biomechanics and biophysics communities.

Numerical analysis and comparison of flow fields in normal larynx and larynx with unilateral vocal fold paralysis

In this study, laryngeal flow fields are investigated and compared in normal larynx and models of larynx with unilateral vocal fold paralysis (UVFP). In paralytic models, three fixed initial glottal gaps are considered to understand the positive or probable negative impacts of surgical operation on unilaterally paralytic larynx, by which the paralyzed vocal fold is brought closer to the mid-plane. Various features of the flow fields have been discussed in detail including glottal gap width, glottal flow rate, glottal exit pressure pattern and glottal jet evolution. The numerical solution of fluid-structure interaction is carried out using ANSYS, and the results confirm some of the favorable effects of surgery on the patient's larynx. It is also shown that by tightening the glottal gap, some of the problems caused by the presence of a motionless vocal fold, such as leakage through glottal gap in the closure phase resulting in breathy voice can be moderated, although some of the symptoms of this disorder remain relatively unchanged.

Vocal fold contact patterns based on normal modes of vibration

The fluid-structure interaction and energy transfer from respiratory airflow to self-sustained vocal fold oscillation continues to be a topic of interest in vocal fold research. Vocal fold vibration is driven by pressures on the vocal fold surface, which are determined by the shape of the glottis and the contact between vocal folds. Characterization of three-dimensional glottal shapes and contact patterns can lead to increased understanding of normal and abnormal physiology of the voice, as well as to development of improved vocal fold models, but a large inventory of shapes has not been directly studied previously. This study aimed to take an initial step toward characterizing vocal fold contact patterns systematically. Vocal fold motion and contact was modeled based on normal mode vibration, as it has been shown that vocal fold vibration can be almost entirely described by only the few lowest order vibrational modes. Symmetric and asymmetric combinations of the four lowest normal modes of vibration were superimposed on left and right vocal fold medial surfaces, for each of three prephonatory glottal configurations, according to a surface wave approach. Contact patterns were generated from the interaction of modal shapes at 16 normalized phases during the vibratory cycle. Eight major contact patterns were identified and characterized by the shape of the flow channel, with the following descriptors assigned: convergent, divergent, convergent-divergent, uniform, split, merged, island, and multichannel. Each of the contact patterns and its variation are described, and future work and applications are discussed.

Computational Models of Laryngeal Aerodynamics: Potentials and Numerical Costs

Human phonation is based on the interaction between tracheal airflow and laryngeal dynamics. This fluid-structure interaction is based on the energy exchange between airflow and vocal folds. Major challenges in analyzing the phonatory process in-vivo are the small dimensions and the poor accessibility of the region of interest. For improved analysis of the phonatory process, numerical simulations of the airflow and the vocal fold dynamics have been suggested. Even though most of the models reproduced the phonatory process fairly well, development of comprehensive larynx models is still a subject of research. In the context of clinical application, physiological accuracy and computational model efficiency are of great interest. In this study, a simple numerical larynx model is introduced that incorporates the laryngeal fluid flow. It is based on a synthetic experimental model with silicone vocal folds. The degree of realism was successively increased in separate computational models and each model was simulated for 10 oscillation cycles. Results show that relevant features of the laryngeal flow field, such as glottal jet deflection, develop even when applying rather simple static models with oscillating flow rates. Including further phonatory components such as vocal fold motion, mucosal wave propagation, and ventricular folds, the simulations show phonatory key features like intraglottal flow separation and increased flow rate in presence of ventricular folds. The simulation time on 100 CPU cores ranged between 25 and 290 hours, currently restricting clinical application of these models. Nevertheless, results show high potential of numerical simulations for better understanding of phonatory process.

Comparison of a fiber-gel finite element model of vocal fold vibration to a transversely isotropic stiffness model

A fiber-gel vocal fold model is compared to a transversely isotropic stiffness model in terms of normal mode vibration. The fiber-gel finite element model (FG-FEM) consists of a series of gel slices, each with a two-dimensional finite element mesh, in a plane transverse to the tissue fibers. The gel slices are coupled with fibers under tension in the anterior-posterior dimension. No vibrational displacement in the fiber-length direction is allowed, resulting in a plane strain state. This is consistent with the assumption of transverse displacement of a simple string, offering a wide range of natural frequencies (well into the kHz region) with variable tension. For low frequencies, the results compare favorably with the natural frequencies of a transversely isotropic elastic stiffness model (TISM) in which the shear modulus in the longitudinal plane is used to approximate the effect of fiber tension. For high frequencies, however, the natural frequencies do not approach the string mode frequencies unless plane strain is imposed on the TISM model. The simplifying assumption of plane strain, as well as the use of analytical closed-form shape functions, allow for substantial savings in computational time, which is important in clinical and exploratory applications of the FG-FEM model.

Fully-coupled aeroelastic simulation with fluid compressibility - For application to vocal fold vibration

In this study, a fully-coupled fluid-structure interaction model is developed for studying dynamic interactions between compressible fluid and aeroelastic structures. The technique is built based on the modified Immersed Finite Element Method (mIFEM), a robust numerical technique to simulate fluid-structure interactions that has capabilities to simulate high Reynolds number flows and handles large density disparities between the fluid and the solid. For accurate assessment of this intricate dynamic process between compressible fluid, such as air and aeroelastic structures, we included in the model the fluid compressibility in an isentropic process and a solid contact model. The accuracy of the compressible fluid solver is verified by examining acoustic wave propagations in a closed and an open duct, respectively. The fully-coupled fluid-structure interaction model is then used to simulate and analyze vocal folds vibrations using compressible air interacting with vocal folds that are represented as layered viscoelastic structures. Using physiological geometric and parametric setup, we are able to obtain a self-sustained vocal fold vibration with a constant inflow pressure. Parametric studies are also performed to study the effects of lung pressure and vocal fold tissue stiffness in vocal folds vibrations. All the case studies produce expected airflow behavior and a sustained vibration, which provide verification and confidence in our future studies of realistic acoustical studies of the phonation process.

Computational Modeling of Fluid-Structure-Acoustics Interaction during Voice Production

The paper presented a three-dimensional, first-principle based fluid-structure-acoustics interaction computer model of voice production, which employed a more realistic human laryngeal and vocal tract geometries. Self-sustained vibrations, important convergent-divergent vibration pattern of the vocal folds, and entrainment of the two dominant vibratory modes were captured. Voice quality-associated parameters including the frequency, open quotient, skewness quotient, and flow rate of the glottal flow waveform were found to be well within the normal physiological ranges. The analogy between the vocal tract and a quarter-wave resonator was demonstrated. The acoustic perturbed flux and pressure inside the glottis were found to be at the same order with their incompressible counterparts, suggesting strong source-filter interactions during voice production. Such high fidelity computational model will be useful for investigating a variety of pathological conditions that involve complex vibrations, such as vocal fold paralysis, vocal nodules, and vocal polyps. The model is also an important step toward a patient-specific surgical planning tool that can serve as a no-risk trial and error platform for different procedures, such as injection of biomaterials and thyroplastic medialization.

The Effect of False Vocal Folds on Laryngeal Flow Resistance in a Tubular Three-dimensional Computational Laryngeal Model

OBJECTIVE

The current study used a three-dimensional (3D) computational laryngeal model to investigate the effect of false vocal folds (FVFs) on laryngeal flow resistance.

METHOD

A 3D, tubular shaped computational laryngeal model was designed with a high level of realism with respect to the human laryngeal anatomy. Two cases, one with the FVFs and the other without the FVFs, were created in the numerical simulation to compare the laryngeal flow behaviors.

RESULTS AND CONCLUSION

The results were discussed in a comparative manner with the previous two-dimensional (2D) computational model. On the one hand, the results demonstrated the similar mechanism as observed in the 2D model that the presence of the FVFs suppressed the deflection of the glottal jet and in doing so, reduced the mixing-related minor loss in the supraglottal region. On the other hand, the 3D flow was more stable and straighter, so the effect of FVFs on suppressing the jet deflection in the 3D model was not as prominent as in the 2D model. Furthermore, the presence of the FVFs also increased the friction-related major loss due to the increased velocity gradient in the restricted flow channel. Therefore, it was hypothesized that the final effect of the FVFs on flow resistance is the combined effect of the reduced mixing-related minor loss and increased friction-related major loss, both of which are highly related to the gap between the FVFs.

Subject-Specific Computational Modeling of Evoked Rabbit Phonation

When developing high-fidelity computational model of vocal fold vibration for voice production of individuals, one would run into typical issues of unknown model parameters and model validation of individual-specific characteristics of phonation. In the current study, the evoked rabbit phonation is adopted to explore some of these issues. In particular, the mechanical properties of the rabbit's vocal fold tissue are unknown for individual subjects. In the model, we couple a 3D vocal fold model that is based on the magnetic resonance (MR) scan of the rabbit larynx and a simple one-dimensional (1D) model for the glottal airflow to perform fast simulations of the vocal fold dynamics. This hybrid three-dimensional (3D)/1D model is then used along with the experimental measurement of each individual subject for determination of the vocal fold properties. The vibration frequency and deformation amplitude from the final model are matched reasonably well for individual subjects. The modeling and validation approaches adopted here could be useful for future development of subject-specific computational models of vocal fold vibration.

Computational modeling of vibration-induced systemic hydration of vocal folds over a range of phonation conditions

Predicting phonation conditions that are benign to voice health remains a biomechanically relevant problem. Our objective is to provide insight into vocal fold (VF) hydration based on continuum-based VF models that are able to compute VF stresses during phonation and a scheme for the extraction and generalization of such computational data based on the principle of linear superposition. Because VF tissue is poroelastic, spatial gradients of VF hydrostatic stresses computed for a given phonation condition determine VF interstitial fluid flow. The present approach transforms, based on linear superposition principles, the computed interstitial fluid velocities at the particular phonation to those at an arbitrary phonation condition. Intersititial fluid flow characteristics for a range of phonation conditions are compared. For phonation conditions with no or moderate collision, no dehydration per vibration cycle is predicted throughout the VF. For more severe collision conditions, tissue dehydration is restricted to a region close to the glottal surface. Interstitial fluid displacement in the VF is found to be heterogeneous and strongly dependent on the phonation condition. A phonation condition is found to exist for which dehydration peaks. The proposed method significantly expands the scope and relevance of conducting isolated numerical simulations of VF vibration.

Validation of a flow-structure-interaction computation model of phonation

Computational models of vocal fold (VF) vibration are becoming increasingly sophisticated, their utility currently transiting from exploratory research to predictive research. However, validation of such models has remained largely qualitative, raising questions over their applicability to interpret clinical situations. In this paper, a computational model with a segregated implementation is detailed. The model is used to predict the fluid-structure interaction (FSI) observed in a physical replica of the VFs when it is excited by airflow. Detailed quantitative comparisons are provided between the computational model and the corresponding experiment. First, the flow model is separately validated in the absence of VF motion. Then, in the presence of flow-induced VF motion, comparisons are made of the flow pressure on the VF walls and of the resulting VF displacements. Self-similarity of spatial distributions of flow pressure and VF displacements is highlighted. The self-similarity leads to normalized pressure and displacement profiles. It is shown that by using linear superposition of average and fluctuation components of normalized computed displacements, it is possible to determine displacements in the physical VF replica over a range of VF vibration conditions. Mechanical stresses in the VF interior are related to the VF displacements, thereby the computational model can also determine VF stresses over a range of phonation conditions.

Feedback control as a framework for understanding tradeoffs in biology

Control theory arose from a need to control synthetic systems. From regulating steam engines to tuning radios to devices capable of autonomous movement, it provided a formal mathematical basis for understanding the role of feedback in the stability (or change) of dynamical systems. It provides a framework for understanding any system with regulation via feedback, including biological ones such as regulatory gene networks, cellular metabolic systems, sensorimotor dynamics of moving animals, and even ecological or evolutionary dynamics of organisms and populations. Here, we focus on four case studies of the sensorimotor dynamics of animals, each of which involves the application of principles from control theory to probe stability and feedback in an organism's response to perturbations. We use examples from aquatic (two behaviors performed by electric fish), terrestrial (following of walls by cockroaches), and aerial environments (flight control by moths) to highlight how one can use control theory to understand the way feedback mechanisms interact with the physical dynamics of animals to determine their stability and response to sensory inputs and perturbations. Each case study is cast as a control problem with sensory input, neural processing, and motor dynamics, the output of which feeds back to the sensory inputs. Collectively, the interaction of these systems in a closed loop determines the behavior of the entire system.

Computational study of effects of tension imbalance on phonation in a three-dimensional tubular larynx model

OBJECTIVES

The present study explores the use of a continuum-based computational model to investigate the effect of left-right tension imbalance on vocal fold (VF) vibrations and glottal aerodynamics, as well as its implication on phonation. The study allows us to gain new insights into the underlying physical mechanism of irregularities induced by VF tension imbalance associated with unilateral cricothyroid muscle paralysis.

METHODS

A three-dimensional simulation of glottal flow and VF dynamics in a tubular laryngeal model with tension imbalance was conducted by using a coupled flow-structure interaction computational model. Tension imbalance was modeled by reducing by 20% the Young's modulus of one of the VFs, while holding VF length constant. Effects of tension imbalance on vibratory characteristic of the VFs and on the time-varying properties of glottal airflow as well as the aerodynamic energy transfer are comprehensively analyzed.

RESULTS AND CONCLUSIONS

The analysis demonstrates that the continuum-based biomechanical model can provide a good description of phonatory dynamics in tension imbalance conditions. It is found that although 20% tension imbalance does not have noticeable effects on the fundamental frequency, it does lead to a larger glottal flow leakage and asymmetric vibrations of the two VFs. A detailed analysis of the energy transfer suggests that the majority of the energy is consumed by the lateral motion of the VFs and the net energy transferred to the softer fold is less than the one transferred to the normal fold.

Subject-specific computational modeling of human phonation

A direct numerical simulation of flow-structure interaction is carried out in a subject-specific larynx model to study human phonation under physiological conditions. The simulation results compare well to the established human data. The resulting glottal flow and waveform are found to be within the normal physiological ranges. The effects of realistic geometry on the vocal fold dynamics and the glottal flow are extensively examined. It is found that the asymmetric anterior-posterior laryngeal configuration produces strong anterior-posterior asymmetries in both vocal fold vibration and glottal flow which has not been captured in the simplified models. It needs to be pointed out that the observations from the current numerical simulation are only valid for the flow conditions investigated. The limitations of the study are also discussed.

Fluid-structure interaction involving large deformations: 3D simulations and applications to biological systems

Three-dimensional fluid-structure interaction (FSI) involving large deformations of flexible bodies is common in biological systems, but accurate and efficient numerical approaches for modeling such systems are still scarce. In this work, we report a successful case of combining an existing immersed-boundary flow solver with a nonlinear finite-element solid-mechanics solver specifically for three-dimensional FSI simulations. This method represents a significant enhancement from the similar methods that are previously available. Based on the Cartesian grid, the viscous incompressible flow solver can handle boundaries of large displacements with simple mesh generation. The solid-mechanics solver has separate subroutines for analyzing general three-dimensional bodies and thin-walled structures composed of frames, membranes, and plates. Both geometric nonlinearity associated with large displacements and material nonlinearity associated with large strains are incorporated in the solver. The FSI is achieved through a strong coupling and partitioned approach. We perform several validation cases, and the results may be used to expand the currently limited database of FSI benchmark study. Finally, we demonstrate the versatility of the present method by applying it to the aerodynamics of elastic wings of insects and the flow-induced vocal fold vibration.

Investigation of prescribed movement in fluid-structure interaction simulation for the human phonation process

In a partitioned approach for computational fluid-structure interaction (FSI) the coupling between fluid and structure causes substantial computational resources. Therefore, a convenient alternative is to reduce the problem to a pure flow simulation with preset movement and applying appropriate boundary conditions. This work investigates the impact of replacing the fully-coupled interface condition with a one-way coupling. To continue to capture structural movement and its effect onto the flow field, prescribed wall movements from separate simulations and/or measurements are used. As an appropriate test case, we apply the different coupling strategies to the human phonation process, which is a highly complex interaction of airflow through the larynx and structural vibration of the vocal folds (VF). We obtain vocal fold vibrations from a fully-coupled simulation and use them as input data for the simplified simulation, i.e. just solving the fluid flow. All computations are performed with our research code CFS++, which is based on the finite element (FE) method. The presented results show that a pure fluid simulation with prescribed structural movement can substitute the fully-coupled approach. However, caution must be used to ensure accurate boundary conditions on the interface, and we found that only a pressure driven flow correctly responds to the physical effects when using specified motion.

The role of finite displacements in vocal fold modeling

Human vocal folds experience flow-induced vibrations during phonation. In previous computational models, the vocal fold dynamics has been treated with linear elasticity theory in which both the strain and the displacement of the tissue are assumed to be infinitesimal (referred to as model I). The effect of the nonlinear strain, or geometric nonlinearity, caused by finite displacements is yet not clear. In this work, a two-dimensional model is used to study the effect of geometric nonlinearity (referred to as model II) on the vocal fold and the airflow. The result shows that even though the deformation is under 1 mm, i.e., less than 10% of the size of the vocal fold, the geometric nonlinear effect is still significant. Specifically, model I underpredicts the gap width, the flow rate, and the impact stress on the medial surfaces as compared to model II. The study further shows that the differences are caused by the contact mechanics and, more importantly, the fluid-structure interaction that magnifies the error from the small-displacement assumption. The results suggest that using the large-displacement formulation in a computational model would be more appropriate for accurate simulations of the vocal fold dynamics.

Influence of numerical model decisions on the flow-induced vibration of a computational vocal fold model

Computational vocal fold models are often used to study the physics of voice production. In this paper the sensitivity of predicted vocal fold flow-induced vibration and resulting airflow patterns to several modeling selections is explored. The location of contact lines used to prevent mesh collapse and assumptions of symmetry were found to influence airflow patterns. However, these variables had relatively little effect on the vibratory response of the vocal fold model itself. Model motion was very sensitive to Poisson's ratio. The importance of these parameter sensitivities in the context of vocal fold modeling is discussed.

Voicing produced by a constant velocity lung source

An investigation is made of the influence of subglottal boundary conditions on the prediction of voiced sounds. It is generally assumed in mathematical models of voicing that vibrations of the vocal folds are maintained by a constant subglottal mean pressure pI, whereas voicing is actually initiated by contraction of the chest cavity until the subglottal pressure becomes large enough to separate the vocal folds. The problem is reformulated to determine voicing characteristics in terms of a prescribed volumetric flow rate Qo of air from the lungs-the evolution of the resulting time-dependent subglottal mean pressure p[overline]-(t) is then governed by glottal mechanics, the aeroacoustics of the vocal tract, and the influence of continued contraction of the lungs. The new problem is analyzed in detail for an idealized mechanical vocal system that permits precise specification of all boundary conditions. Predictions of the glottal volume velocity pulse shape are found to be in good general agreement with the traditional constant-pI theory when pI is set equal to the time averaged value of p[overline]-(t). But, in all cases examined the constant-pI approximation yields values of the mean flow rates Qo and sound pressure levels that are smaller by as much as 10%.