Elsevier

Journal of Computational Physics

Volume 227, Issue 22, 20 November 2008, Pages 9303-9332
Journal of Computational Physics

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

https://doi.org/10.1016/j.jcp.2008.05.001Get rights and content

Abstract

A new numerical approach for modeling a class of flow–structure interaction problems typically encountered in biological systems is presented. In this approach, a previously developed, sharp-interface, immersed-boundary method for incompressible flows is used to model the fluid flow and a new, sharp-interface Cartesian grid, immersed-boundary method is devised to solve the equations of linear viscoelasticity that governs the solid. The two solvers are coupled to model flow–structure interaction. This coupled solver has the advantage of simple grid generation and efficient computation on simple, single-block structured grids. The accuracy of the solid-mechanics solver is examined by applying it to a canonical problem. The solution methodology is then applied to the problem of laryngeal aerodynamics and vocal fold vibration during human phonation. This includes a three-dimensional eigen analysis for a multi-layered vocal fold prototype as well as two-dimensional, flow-induced vocal fold vibration in a modeled larynx. Several salient features of the aerodynamics as well as vocal fold dynamics are presented.

Introduction

Flow–structure interaction (FSI) is a common phenomenon in biological systems. Typical examples related to biomedical engineering include the cardiovascular system (heart valves and arteries), and the larynx. The ability to computationally model the flow–structure interaction in these systems could help us understand the underlying biophysics, investigate pathologies, and potentially advance medical treatments. Structural flexibility and flow–induced deformation is also ubiquitous in nature. For instance, flow–structure interaction is a key feature in biological locomotion including fish/mammalian swimming [1] and insect/bird flight, and the ability to model this interaction is important in learning the underlying physics of these modes of locomotion.

One of the main challenges in developing such biophysical models is handling of the complex and moving anatomical geometries. The finite-element method (FEM) is the traditional way of dealing with complicated computational domains (e.g. [2], [3]). However, grid generation and solution of the associated algebraic equations can be quite expensive. Furthermore, biological configurations present a singularly difficult proposition for such methods given the highly complex geometries, motions, deformation and material properties that are usually encountered in these configurations.

In recent years, the immersed-boundary (IB) method has gained popularity in computational fluid dynamics (CFD) for handling complex and/or moving boundaries. In the IB method, a structured, usually Cartesian, grid which does not conform to the flow boundary is used for discretizing the governing equations [4]. Recent review on the IB method and its variants can be found in Mittal and Iaccarino [5]. Compared to the boundary-conforming structured and unstructured methods, the IB method has the advantages of simple grid generation [4] and ease of incorporating multigrid [6] and domain-decomposition based parallel algorithms [7].

The Cartesian grid based IB method has also been applied in the computation of solid-mechanics. For example, Sethian and Wiegmann [8] used a type of IB method to solve linear elastostatics on arbitrary two-dimensional domains and the solution was used in an optimization procedure to iteratively improve structural design. In their approach, a level-set method was used to represent the boundaries of the solid body, and an immersed-boundary method based on Li and LeVeque [9] and Li [10] was used to prescribe the discontinuities in the governing equations across the solid/void boundary. This approach allowed them to change the geometry and topology of the structure during the optimization process without modifying the underlying grid.

Udaykumar and coworkers [11], [12] used an Eulerian method to simulate high-speed multi-material impact. Their method was based on a fixed Cartesian grid and a sharp-interface IB method was used to deal with large deformations of the material–material and material–void interfaces. The approach was particularly attractive in that the issues associated with severe mesh distortion and entangling, which would be faced by conventional body-conformal methods, can be circumvented.

In this paper, we present a Cartesian grid based approach for modeling a class of FSI problems typically encountered in biological applications. More specifically, we employ the previous sharp-interface IB method [13], [7], [14] to solve the Navier–Stokes equations that govern the flow, and devise a new IB formulation that allows us to compute the linear elastodynamics of complex three-dimensional (3D) structures. FSI is accomplished by operating the two solvers in a coupled manner. Compared to the IB methods described in [8], [11], [12], our method can be used for simulating dynamics of linearly elastic or viscoelastic solids as well as flow-induced deformation of such solids. The FSI solver is also designed to solve two- as well as three-dimensional problems and is therefore very well suited for high-fidelity modeling of biological configurations.

Although the IB method we present here for the 3D linear viscoelasticity is inspired from the approach developed in the context of the fluid dynamics by Mittal and coworkers [13], [7], [14] and therefore bears some similarity to that approach, the implementation is significantly different, especially with regard to the treatment of the traction boundary condition which is a unique feature of solid dynamics. This issue is discussed in detail in Section 2. It should also be noted that this method is different from the IB method described in Li and coworkers [9], [10] and Sethian and Wiegmann [8]. In their methods, the solution experiences discontinuities across the singular interface immersed in the domain, and the finite-difference formulas involving the nodes across the interface were corrected by using Taylor’s series around the interface and taking into consideration of the discontinuities. In contrast, our method is based on a ghost-cell methodology where the ghost-node value is a smooth extrapolation from the solution on the physical side of the boundary. There is no discontinuity involved at the boundary in our method. Furthermore, those methods require derivation of the correction term in the finite-difference formulas near the boundary, which in our view is inconvenient if applied to the 3D elasticity. In comparison, the finite-difference equations in our method are standard formulations and are thus much simpler.

Finally, the current method differs from the extended IB method or immersed finite-element method proposed in [15], [16] in that, in our formulation, (1) there is no body force imposed at the fluid/solid boundary or within the solid body, (2) only Cartesian meshes are used.

A particular focus of the current work is developing a computational modeling capability that can capture the physics of phonation which refers to the process of sound production in the larynx. Phonation is essentially a result of flow-induced vibration of the vocal folds (VF). Fig. 1(a) shows a coronal (front-to-back) view of larynx obtained from a computed tomography (CT) scan. The image clearly shows the two vocal folds that protrude into the airway inside the larynx. During phonation, the two VFs are brought together at the midline and tightened so as to obstruct the passage of air from the lungs to the vocal tract above. Air is then forced through this laryngeal passageway (called the “glottis”) due to buildup of the pressure inside the lungs, and this results in sustained flow-induced vibration of the VFs during which air is expelled into the vocal tract as an oscillatory jet called the “glottal jet”. This sustained vibration and the oscillating airflow give rise to the generation and propagation of sound, and this process is called phonation. When attention is focused on the VF vibration and the jet behavior, air compressibility is often neglected and an incompressible flow can be assumed. The vocal fold itself has a complex structure as shown in Fig. 1(b), and the various constituents of the VF are known to play distinct roles in the vibratory dynamics [18].

The dynamics of the vocal folds and glottal jet are difficult to examine in experiments. Thus, despite a significant number of in vitro and in vivo studies [18], much remains to be understood regarding the biophysics of phonation. For example, little has been known about the unsteady vortex motions in the supraglottal region and their effect on the vocal fold vibration as well as the sound generation. A mathematical model that describes the dynamical process of phonation could complement experimental studies thereby helping us understand the physics of voice production. It may also have potential significance for examining certain voice pathologies and treating voice disorders.

In the past, a number of mathematical models of different complexity have been developed for describing the FSI process of phonation. The first study that attempted to examine the phonation physics using a fully coupled FSI approach was that of Ishizaka and Flanagan [19]. In this pioneering study, the VFs were modeled as two lumped masses and the air was treated as a one-dimensional inviscid flow. The model was extremely simple but was able to successfully demonstrate sustained flow-induced vibrations. Subsequent to this, lumped-mass models with more degrees of freedom were proposed and employed, e.g., the sixteen-mass model used by Titze [20]. Low-order models have been used with varying degrees of success to study some specific features of phonation. For example, the chaotic behavior in the VF vibration was examined by Jiang et al. [21] using a two-mass model.

Continuum models of the vocal folds have been employed in recent years. A two-/three-dimensional hybrid FEM model of VFs was introduced by Alipour et al. [22] where the VF tissues were assumed to have three layers and each layer was transversely isotropic and governed by linear viscoelasticity. Coupling this model with a two-dimensional (2D) flow solver, Alipour and Scherer [23] studied the bulging effect of the medial surface of the VFs due to glottal adduction. Rosa et al. [24] presented a fully 3D model in which dynamics of the three-layer and transversely isotropic VFs was coupled with an incompressible flow solver to simulate the FSI. In addition, they included the contact force during the VF closure, and the false vocal folds and laryngeal ventricles were incorporated into their simulation to better approximate the physical geometry. Both the solid dynamics and the fluid dynamics were solved using the FEM methods. Using the model, the authors examined the phase difference in the VF tissue deformation and the effect of the false VFs on the pressure distribution over the laryngeal surfaces.

Tao and Jiang [25] also recently considered a 3D VF model. Combining the model with Bernoulli’s law, they investigated the anterior–posterior biphonation (simultaneous occurrence of two independent fundamental frequencies during phonation) phenomenon. Thomson et al. [26] used both 2D FEM simulations and experiments on a synthetic VF model to study the energy transfer from the airflow to the VF during the FSI. Hunter et al. [27], [28] used numerical simulations to describe the dynamics of the VF abduction and adduction – the posturing movements of the VFs during phonation aside from their vibration.

All of the above models have been useful in describing the basic vibratory function of the VFs and some particular aspects of the phonation process. However, for a more detailed analysis of phonation, higher-fidelity models that can incorporate more realistic geometries and provide higher accuracy both in the fluid and solid dynamics are needed. Furthermore, in order to examine patient-specific configurations, which is key to effective treatment, an efficient method is needed for rapid modeling of a variety of configurations. In our research, we attempt to develop a continuum mechanics based methodology which can resolve a large range of temporal and spatial scales in both the VF vibration and aerodynamics. This model will be able to capture details of the vibratory characteristics as well as the flow behavior, and thus allow us to gain a deeper insight into the physics of phonation. The model is expected to eventually be used for improving the outcome of laryngeal surgeries. For example, in medialization laryngoplasty, a surgical procedure used to treat vocal fold paresis and paralysis, a uniquely configured structural implant is inserted into the diseased VF to improve its vibratory characteristics [29]. A high-fidelity computational model could potentially help surgeons predict the effect of the implant and possibly improve the success rate of this procedure [30]. This indeed is the long-term goal of the current effort.

In this paper we describe a crucial step toward that goal. We have developed a new Cartesian grid based immersed-boundary method to simulate the elastodynamics of complex elastic and viscoelastic solid structures. This solver is coupled with an existing IB method that solves the incompressible Navier–Stokes equations. This combined method allows us to model FSI with complex geometries with relative ease. In Section 2, we describe the IB method for general viscoelastic solids subject to linear deformation. The method is validated and its accuracy is tested using a canonical problem and the grid refinement in Section 3. In Section 4.1, we apply the IB method to the problem of phonation and compute the vibration modes of a prototypical 3D VF. In Section 4.2, we couple the method with an immersed-boundary flow solver to simulate the flow-induced VF vibration in two dimensions. Summary and conclusions are given in Section 5.

Section snippets

An immersed-boundary method for linear viscoelasticity

In the following, we describe the salient features of the numerical method developed to solve the dynamical equations of a linear viscoelastic solid. We first describe the underlying methodology for solving the governing equations on a Cartesian mesh and then describe how the appropriate boundary conditions are applied over the immersed boundaries that do not conform to the Cartesian mesh.

Grid refinement study

The spatial accuracy of the linear-elastic solver as well as its fidelity is examined by computing the numerical solution for a non-trivial geometry on different grids and comparing with a known exact solution. Here we consider an infinitely long annulus with inner radius R1 and outer radius R2 as shown in Fig. 6(a). The outer surface of the annulus is displaced in the radial direction by distance s, and the inner surface is either fixed (i.e., zero displacement) or free (i.e., zero traction).

Application to Phonation

Histologically, the vocal fold consists of the vocalis muscles and mucosa, and the mucosa is comprised of the epithelium at the surface and the lamina propria below, as shown in Fig. 1(b). At the VF edge, the lamina propria can be further divided into three layers: the superficial, intermediate, and deep layers. From a mechanical point of view, these layers may be regrouped into three layers: the cover (the epithelium and superficial layer of the lamina propria); the ligament (the intermediate

Conclusions

We have developed a numerical approach to simulate a class of flow–structure interaction typically encountered in biological systems. A new sharp-interface IB methodology has been developed to solve the elastodynamics of a linear viscoelastic solid. The key feature in the development of the IB methodology for elastodynamics is a robust and efficient formulation for imposing the traction boundary condition on the solid surface. This newly developed elastodynamic solver is coupled to an existing

Acknowledgment

We want to acknowledge Dr. Haibo Dong for his help in understanding the flow solver. This research is supported by NIDCD grant R01 DC007125-01A1. Support from the GW Institute of Biomedical Engineering is also gratefully acknowledged.

References (54)

  • F. Alipour et al.

    Vocal fold bulging effects on phonation using a biophysical computer model

    J. Voice

    (2000)
  • T.D. Anderson et al.

    Thyroplasty revisions: frequency and predictive factors

    J. Voice

    (2003)
  • T. Ye et al.

    An accurate cartesian grid method for simulation of viscous incompressible flows with complex immersed boundaries

    J. Comput. Phys.

    (1999)
  • S. Marella et al.

    Sharp interface cartesian grid method I: an easily implemented technique for 3D moving boundary computations

    J. Comput. Phys.

    (2005)
  • F. Alipour et al.

    Numerical simulation of laryngeal flow in a forced-oscillation glottal model

    Comput. Speech Lang.

    (1996)
  • F. Alipour et al.

    Dynamic glottal pressures in an excised hemilarynx model

    J. Voice

    (2000)
  • I.R. Titze

    Mechanical stress in phonation

    J. Voice

    (1994)
  • R. Mittal

    Computational modeling in biohydrodynamics: trends, challenges, and recent advances

    IEEE J. Oceanic Eng.

    (2004)
  • M. Bathe et al.

    A fluid-structure interaction finite element analysis of pulsatile blood flow through a compliant stenotic artery

    J. Biomech. Eng.

    (1999)
  • R. Mittal et al.

    Immersed boundary methods

    Annu. Rev. Fluid Mech.

    (2005)
  • R.J. LeVeque et al.

    The immersed interface method for elliptic equations with discontinuous coefficients and singular sources

    SIAM J. Numer. Anal.

    (1994)
  • Z. Li

    A fast interactive algorithm for elliptic interface problems

    SIAM J. Numer. Anal.

    (1998)
  • H. Dong et al.

    Wake topology and hydrodynamic performance of low-aspect-ratio flapping foils

    J. Fluid Mech.

    (2006)
  • M. Hirano et al.

    Histological Color Atlas of the Human larynx

    (1993)
  • I.R. Titze

    Principles of Voice Production

    (1994)
  • K. Ishizaka et al.

    Synthesis of voiced sounds from a two-mass model of the vocal cords

    Bell Syst. Tech. J.

    (1972)
  • I.R. Titze

    The human vocal cords: a mathematical model, Part I

    Phonetica

    (1973)
  • Cited by (172)

    View all citing articles on Scopus
    View full text