Journal of Biomedical Science and Engineering
Vol.7 No.3(2014), Article ID:43331,16 pages DOI:10.4236/jbise.2014.73018

Modeling of Soft Tissues Interacting with Fluid (Blood or Air) Using the Immersed Finite Element Method

Lucy T. Zhang

JEC 2049, Department of Mechanical, Aerospace, and Nuclear Engineering, Troy, USA


Copyright © 2014 Lucy T. Zhang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. In accordance of the Creative Commons Attribution License all Copyrights © 2014 are reserved for SCIRP and the owner of the intellectual property Lucy T. Zhang. All Copyright © 2014 are guarded by law and by SCIRP as a guardian.

Received 26 November 2013; revised 30 December 2013; accepted 7 January 2014

KEYWORDS:Biomechanics; Fluid-Structure Interactions; Biomechanics


This paper presents some biomedical applications that involve fluid-structure interactions which are simulated using the Immersed Finite Element Method (IFEM). Here, we first review the original and enhanced IFEM methods that are suitable to model incompressible or compressible fluid that can have densities that are significantly lower than the solid, such as air. Then, three biomedical applications are studied using the IFEM. Each of the applications may require a specific set of IFEM formulation for its respective numerical stability and accuracy due to the disparities between the fluid and the solid. We show that these biomedical applications require a fully-coupled and stable numerical technique in order to produce meaningful results. Biomechanics; Fluid-Structure Interactions; Biomechanics


In the past decade, the interest in developing novel simulation techniques for modeling fluid-structure interactions revived due to the increasing demands in capabilities to accurately and efficiently study biomedical applications. Biomedical applications often involve fluid (blood) interacting with soft tissues. In some cases, the fluid can also be air, which has a disparate density compared to the soft tissues.

Since soft tissues come in with all forms, shapes and sizes, it is more convenient to set up a simulation using non-boundary-fitted modeling technique. The non-boundary-fitted approaches avoid the re-meshing process by defining independent meshes for the fluid and solid respectively. The solid can freely move on top of the fluid grid without deforming the surrounding fluid. A widely used numerical approach for bio-interface applications is the immersed boundary (IB) method, which was initially proposed by Peskin to study the blood flow around heart valves [1-7]. The mathematical formulation of the IB method employs a mixture of Eulerian and Lagrangian descriptions for fluid and solid domains. In particular, the entire fluid domain is represented by a uniform background grid, which can be solved by finite difference method with periodic boundary conditions; whereas the submerged structure is represented by a fiber or boundary network. The interaction between the fluid and structure is accomplished by distributing the nodal forces and interpolating the velocities between Eulerian and Lagrangian domains through a smoothed approximation of the Dirac delta function. The advantage of the IB method is that the fluid-structure interface is tracked automatically by following the displaced structural boundary movement, which removes the costly computations due to various mesh update algorithms. Many other numerical algorithms have been developed that are inspired by the IB method, such as the immersed interface method (IIM) [8-14], the extended immersed boundary method (EIBM) [15] and the immersed boundary finite element method (IB-FEM) [16,17]. A review on several methods can be found in [18].

The problem existing in the non-boundary-fitted approaches mentioned above is the lack of more realistic representations of the solid, which hinders the accurate assessment of the material behavior and its deformation. Since the solid and the fluid domains are fully-coupled, a slight inaccuracy in estimating the solid solution may even affect the surrounding fluid solutions. This problem may propagate over time and cause instabilities in the solution or convergence issues. The immersed finite element method (IFEM) [19-22] was developed to tackle this problem by representing the background viscous fluid with an unstructured finite element mesh and non-linear finite elements for the immersed deformable solid. Similar to the immersed boundary method, the fluid domain is defined on a fixed Eulerian grid. However, the solid domain is constructed independently with a Lagrangian mesh, which makes it possible to use a more detailed constitutive model to describe the solid material such as linear elastic, hyperelastic and viscoelastic. This approach is particularly attractive to modeling biomedical applications including stent deployment, blood flow in athersclerosis in arteries, etc. [23-31]. With finite element formulations for both fluid and solid domains, the submerged structure is solved more realistically and accurately in comparison with the corresponding fiber network representation in the IB method. The caveat, of course, is to remove the artificial fluid where the solid volume occupies. Since the solid moves at every time step, the artificial fluid also moves. Using the non-boundary-fitted approach, this volume can be easily identified. The fluid solver is based on a stabilized equal-order finite element formulation applicable to problems involving moving boundaries [32-34]. This stabilized formulation prevents numerical oscillations without introducing excessive numerical dissipations. It is also possible to assign sufficiently refined fluid mesh in local regions wherever necessary to obtain more accurate interfacial solutions.

The two-way coupled approach, i.e. the interpolation and the distribution of the velocity and the forces between the two domains, is quite robust when the solid behaves very much like the fluid. However, if there exists high discontinuity in density as well as other intrinsic parameters in the solid and the fluid, the force and the velocity to be interpolated between the two fields can no longer provide consistent convergence. Therefore, a semi-implicit algorithm for the immersed finite element [35] was developed to alleviate the situations when large density difference and/or stiff solid material are used in the solid domain. The calculation of the fluid-structure interaction force is modified in order to achieve a larger stabilization range.

The semi-implicit IFEM algorithm works well when the fluid dominates the dynamics of the system, in which the solid moves and deforms by following the fluid flow. However, when the solid dynamics or inertia must be taken into account, letting the solid follow the fluid movement may lead to unrealistic solid deformation and sometimes even causes the distortion of the solid mesh, because it is not appropriate to still approximate the solid behaviors using only the fluid velocity. A modified IFEM algorithm (mIFEM) [36] was then introduced to provide a more accurate prediction of the solid motion and deformation, which directly depend on the solid inertia effects, constitutive laws, and the fluid solutions near the fluidstructure interface. The results show that it produces more accurate and reasonable solid responses compared to the original IFEM algorithm.

In this paper, we first provide a detailed derivation and descriptions of the mathematical formulations for IFEM, the semi-implicit IFEM, and the modified IFEM. Then, we will show several biomedical applications that used IFEM algorithms. The choice of the algorithm used for each application is dependent on the nature of the fluidstructure interactions involved, which comes clear as the readers go through the motivations of each of the algorithm development.


2.1. Kinematics

Let us consider a deformable structure that occupies a finite domain, , which is completely immersed in a fluid domain, as illustrated in Figure 1. The fluid and the solid together occupy the entire computational domain Ω, and they intersect at a common interface, where “FSI” represents a line if Ω is a two-dimensional domain or a surface if Ω is in three-dimensions. The interface coincides with the solid boundary. The nomenclature involved can be partitioned into two categories: one belongs to the solid and the other to the fluid. The notations associated with the solid have superscript s to distinguish them from those of the fluid f.

2.2. Assumptions

Before showing derivations, we first need to state three assumptions:

1) The fluid exists everywhere in the domain, Ω. This assumption allows us to generate fluid and solid meshes and solve fluid and solid equations independently, thus avoiding frequent mesh updating schemes required to track the fluid-structure interface. In the IFEM, the solid immersed in the fluid domain occupies a physical space or volume in the computational domain. Therefore, when the solid domain, , is constructed, it overlaps with the entire domain Ω filled with fluid. Since both the solid and the “artificial” fluid co-exist in, it is also referred to as the “overlapping domain”, , i.e., , in later text. This is illustrated in Figure 1. This assumption may simplify the computations, but does not comply with the actual physics. Therefore, this “artificial” fluid effect in the solid domain must be eliminated when formulating the equations.

2) The interface between the fluid and the solid must abide by the matching velocity (or no-slip) and traction boundary conditions. This assumption states that the sol-

Figure 1. Computational domain decomposition.

id boundary moves together with the artificial fluid boundary or vice-versa, and the surface traction on both domains are equal and opposite. This assumption allows appropriate coupling to occur between the two domains.

3) The solid must always remain immersed in the fluid domain to avoid inaccurate interpolations at the fluidsolid interface.

2.3. Interpolations between Fluid and Solid Domains

The solid and fluid meshes are constructed independently, therefore it is impossible to have the moving solid boundary nodes exactly coinciding with the fluid nodes in Ω. An interpolation function, , must be used to couple the fluid velocity field and the solid nodal velocity, such that:


Similarly, the interaction force calculated in the solid (overlapping) domain is distributed to the fluid domain as:


This two-way coupling is necessary to ensure the stability and convergence of the algorithm. There are multiple ways of performing the interpolation process. The interpolation function, , can be acquired through the discretized Dirac delta function [37], the sharp finite element interpolation function [22], or the reproducing kernel interpolation function [20]. The details and the characteristics of each approach can be found and compared in Ref. [22].


The derivation of IFEM starts from the principle of virtual work or the weak form, which is used for standard finite element analysis. The weak forms of the derived equations are equivalent to their strong forms if the weak form solution is smooth enough to satisfy at least continuity.

3.1. Derivation

The virtue work in the solid domain with a test function, can be expressed as:


The terms in the bracket of Equation (3) describe the governing equation for the solid, where is the stress which is directly related to the internal force that is determined by the material types and properties. The term, or, is the inertial force and is the body or external force.

To include the artificial fluid related terms without contradicting the equilibrium, Equation (3) can be rewritten as:


The added terms are underlined. One can notice that they sum to zero. Since and belong to the same physical space, we can rearrange this equation to yield:


Equation (5) now contains two integral terms. The first term is the work done by the solid in the solid domain subtracting the work done by the artificial fluid. The second term represents the work done by the artificial fluid in this overlapping domain.

We define the first term in Equation (5) to be the interaction force:


Since this interaction force is first evaluated in the solid domain on the solid nodal points, it is therefore labeled as. It represents the interaction force acting on the solid from the fluid. Once it is evaluated, the nodal forces are then distributed onto the fluid domain as. This force, then, becomes the driving force for the fluid. Equation (5) becomes:


Using the no-slip assumption made in Assumption (2) which allows in, Equation (7) becomes:


Inside the bracket of this equation is the momentum equation for the artificial fluid, where all the terms resembles the Navier-Stokes equation except the addition of the interaction force. The interaction force only exists in the overlapping region and its immediate surroundings. Its value diminishes to zero at places outside the region.

Now, combining the work done by the real fluid and the artificial fluid described in Equation (8) with the expansion of the total time derivative term, , we obtain


The two integral terms in Equation (9) can be combined into the entire computational domain, , as:


Since the fluid is homogenous and both physical and artificial fluids are assumed to be incompressible, we can write the complete governing equations of fluid as



In the IFEM formulation, the term can be interpreted as the external force applied to the fluid that is generated from the artificial fluid. It is important to note that since the solid nodal velocities follow that of the overlapping fluid grid velocities, the compressibility of the solid must follow that of the fluid as well. Therefore, the solid must be incompressible or at least nearly incompressible when the fluid is incompressible. This restriction is alleviated in the modified IFEM algorithm in Section 5.

3.2. Outline of the IFEM Algorithm

An outline of the IFEM algorithm can be illustrated as follows:

1) Given the structural configuration and the fluid velocity from the previous time step2) Evaluate the nodal interaction forces on solid material points, using Equation (6)3) Distribute the material nodal force onto the fluid grid, from to using interpolation function Equation (2)4) Solve for fluid velocities and pressure implicitly using Equations (11) and (12) at current time step5) Interpolate the velocities in the fluid domain to the material points, i.e. from to, as in Equation (1)6) Update the positions of the structure using and go back to step (1).


In the IFEM, small time step has to be used to ensure the stability of the coupling procedure because the solid domain and fluid domain are coupled to each other explicitly at every time step. Since the Navier-Stokes equations are solved implicitly, such small time step requirement due to the coupling stability makes the whole algorithm numerically inefficient especially for the cases when the solid properties are very different from the fluid. Semi-implicit coupling between the fluid and solid domain is then introduced in order to enlarge the stability range.

4.1. Explicit Fluid-Structure Interaction Force

Although the force or the work is balanced seamlessly in the strong and weak forms at every time step, the fluid domain is numerically balanced with the fluid-structure interaction force evaluated based on the solid configuration of the previous time step. Therefore, the coupling between the two domains is considered explicit. In Equation (6) both the acceleration term and the solid internal stress term are evaluated based on the solid nodal velocity which is interpolated from the fluid velocity of the previous time step,:


with two of the terms evaluated from time step, the interaction force is effectively. It is then passed onto the Navier-Stokes momentum equation to solve for and at the current time step:


Noting that every term in Equation (14) is solved in the current time step except the last term where the interaction force is evaluated from the previous time step. The last term can be related to the current interaction force by taking the Taylor's expansion at time step n:


The error due to the explicit coupling can be approximated by substituting Equation (15) into Equation (14):


Based on the definition of the fluid-structure interaction force in Equation (13), each term of this fluidstructure interaction force contributes to the accumulative coupling error. These terms are proportional to the density ratio, the stiffness ratio and the gravity ratio, respectively. Here, K is an equivalent Young’s modulus of the solid representing the stiffness of the solid material. If any of these terms are large, the resulting error due to the coupling would be large. These large errors often result in instability or divergence of the solution.

4.2. Semi-Implicit Fluid-Structure Interaction Force

To alleviate the numerical issues caused by the restrictions in time step size of explicit coupling and the convergence problem due to highly disparate properties between the fluid and the solid domains, a semi-implicit approach is introduced [35]. In the semi-implicit algorithm, the interaction force is re-defined in the solid domain, which only includes the internal forces for the fluid and solid from the original definition in Equation (13), such that:


The rest of the terms in the original explicit formulation Equation (13), namely, the inertial and external force terms, are now incorporated into the fluid equations. The newly defined is distributed to the fluid domain as in the original IFEM. The Navier-Stokes equations now must also be re-defined as follows:



where is defined as:


Here, the indicator function, , is to identify the real fluid region, the artificial fluid region or the solid region, and the fluid-structure interface, in the computational domain. The value of the indicator function is ranged between 0 and 1 where it is 0 if an entire element belongs to the fluid and 1 if an entire element belongs to the solid. This newly revised fluid’s momentum equation combines the inertial and the gravity terms in the original FSI force equation.

To further improve the algorithm, we enhance the indicator function so that it can accommodate high density ratios between the fluid and the solid. This means that the interfacial elements have the transitional indicator values. To do so, the elements that contain the fluid-solid interface have varying indicator value that transits from 0 to 1 by solving Poisson’s equation [38]:


where is interpolated by


Here, is the unit outward normal of the solid interface and is the same interpolation function used by the velocity interpolation function and force distribution. The boundary conditions are given as,



Since the fluid-solid interface moves, this indicator function is updated at every time step based on the relative position of the solid domain in the entire computational domain. Comparing to the original IFEM algorithm, the inertial and the external force terms in the original interaction force formulation are now been considered in the governing equation and can be evaluated iteratively with the most updated velocity field. This is, therefore, considered as semi-implicit IFEM algorithm.

In this formulation, the solid internal stress is still dependent on the fluid solutions from the previous time step. Even though this term is evaluated explicitly, the semi-implicit scheme still significantly improves the convergence of the solution when the solid material is very stiff. If we re-visit the coupling error Equation (16), the magnitude of the coupling error in the internal stress term is proportional to the stiffness ratio. Noting that in the semi-implicit form, is defined in Equation (20); the stiffness ratio here is in fact. For most of the cases, the solid density will be larger than the fluid density, which reduces the coupling error comparing to from the original explicit form. Therefore, although the solid internal force is still computed explicitly, the coupling error for the semi-implicit scheme is smaller than the explicit scheme.

Overall, this semi-implicit scheme relaxes the small time step requirement and ensures the stability of the FSI force estimation. In particular, this algorithm can handle a much larger range of fluid and solid properties without sacrificing the computational time.

4.3. Outline of the Semi-Implicit IFEM Algorithm

An outline of the semi-implicit IFEM algorithm is illustrated as follows:

1) Given the structural configuration and the fluid velocity from time step2) Evaluate the nodal semi-implicit interaction forces on the solid material points, using Equation (17)3) Distribute the material nodal force onto the fluid grid, from to using delta function Equation (2)4) Obtain the indicator field by solving Equation (21)5) Solve for fluid velocities and pressure implicitly using Equations (18), (19) and (20)6) Interpolate the velocities in the fluid domain to the material points, i.e. from to, as in Equation (1)7) Update the positions of the solid nodes using and go back to step (1).


Although the semi-implicit coupling scheme alleviates some convergence issues when the fluid and solid have large property differences, one can notice that the solid dynamics is still been controlled by the artificial fluid. For cases where the solid behavior dominates the entire system or high Reynolds number flows using the IFEM algorithm may lead to unrealistic solid deformation and may even cause the severe distortion of the solid mesh, because it is not appropriate to approximate the solid behavior based on the fluid velocity by letting. The idea of the modified IFEM is to let the artificial fluid to behave more like the solid, or letting. Doing so allows the solid governing equation to be solved rather than be evaluated. Since the artificial fluid is not real anyway, its role is to produce the same velocity as the solid so that the real fluid realizes the existence of the solid. This modified IFEM algorithm allows the solid behaviors to be estimated more accurately and have stronger influences in the fluid-structure interactions. The detailed rationale and derivations, as well as validation cases were presented in [36].

5.1. Derivation

In order to find the solid displacement field and the velocity field, the solid equation is solved,


The solid stress is evaluated using the solid strain tensor,


where. Different combinations of and provide various choices of solid material constitutive laws such as linear elastic, viscolinear elastic, hyper-elastic, etc.

The boundary condition of the solid domain can be applied using either Dirichlet boundary condition described in Equation (26) or Neumann boundary condition described in Equation (27).



Here, n is the outward normal of the fluid-structure interface. These boundary conditions are evaluated based on the fluid velocity and stress on the fluid-structure interface solved from the fluid equations at previous time step. is the time step size.

Once the solid solution is obtained, the next step is to make the artificial fluid to follow the solid, i.e. solving the artificial fluid governing equation so that in. To accomplish this, the artificial fluid property, such as the density, should mimic that of the solid.

The continuity equation of the artificial fluid can be written as,


It is also necessary that the artificial fluid has the same compressibility as the solid. Therefore, the artificial fluid is described as a pseudo-compressible fluid:


where the compressibility of the solid is used for the compressibility of the artificial fluid. The continuity equation in the artificial fluid domain can be eventually written as follows,


Using the same semi-implicit interaction force definition and the indicator function as mentioned in the semiimplicit IFEM algorithm, the continuity and momentum equations of the fluid domain, which combines the real fluid domain and artificial fluid domain, can be written as follows,



To enforce the assumption in, a correction force is introduced and added into the fluid-structure interaction force. The correction force, is defined as,


The correction force is effectively the difference between the material derivative of velocity in the solid and the artificial fluid so that both the inertial and convective acceleration forces are accounted for. It would be zero if the artificial fluid follows the solid exactly. Including this correction force the fluid structure interaction force is redefined as,


5.2. Outline of the Modified IFEM Algorithm

The algorithm of the modified IFEM is outlined as the following:

1) Solve the solid governing equation Equation (24) with the boundary conditions interpolated from the fluid field in the previous time step2) Evaluate the nodal semi-implicit interaction forces on the solid material points, using Equation (34)3) Distribute the material nodal force onto the fluid grid, from to using interpolation function, Equation (2)4) Obtain the indicator field by solving Equation (21)5) Solve for fluid velocities and pressure implicitly using Equations (31) and (32)6) Interpolate the interface velocities and stress from the fluid domain to the material points, Equations (26) and (27), go back to step (1).


In this paper, three biomedical applications are demonstrated. The first example is a blood cell traveling through a bifurcated blood vessel; the second example is to simulate the deployment of an angioplasty stent, which was first presented in [31]; the third example is to study the vocal folds vibration. The first two examples used the original IFEM algorithm where the fluid is blood and the solid is soft tissues. We used the mIFEM algorithm for the third example where the density ratio between the fluid (air) and the soft tissue is large.

6.1. Red Blood Cells in a Bifurcated Vessel

Understanding the behavior of red blood cell (RBC) flowing in blood vessels, especially when bifurcation happens, is important in estimating the nonuniform hematocrit distribution that would affect the microvasular oxygen distribution, the effective viscosity of blood in microvessels and the distribution of other metabolites. Learning the behaviors of diseased RBCs that have abnormal rigidity, radius and shape, can be helpful in designing medical therapy. Using the established IFEM method, we can simulate the motion and the deformation of the red blood cell within the vessels, and study in detail how the geometry of bifurcation and fluid field affect and direct which daughter branch the RBC flows.

The geometry of a bifurcated blood vessel is shown in Figure 2, where is the diameter of the mother vessel; and are the diameters of the daughter vessels on the top and the bottom, respectively; and are the respective branching angles of the daughter vessels; the branching fillet radii, and are given as to make the vessel branching transition smooth;, and represent the flow rate of each vessel. A RBC is placed near the inlet of the vessel. The radius of the RBC is also given as 2.66 μm. The incoming velocity of the mother vessel is set to be a constant as 0.1 cm/s. The branching angles and are set to be equal and constant,. In this study, we set the diameter of the mother branch to be, and consider two sets of diameter ratios, , to be 1 and 1.44. When the diameter ratio is 1, it is considered as symmetric bifurcation; when it is not 1, then it is considered as asymmetric.

Figures 3 and 4 represent the blood cell behaviors when encountering bifurcation in symmetric and asymmetric vessels with different original positions and ratios of flow rate of each daughter vessel. Based on Figure 3 we can notice that although the daughter branches are symmetric in the geometry, due to the asymmetric boundary conditions where the ratio of the flow rates for the two daughter branches is, the blood cell tends to move to the daughter branch with a higher mass flow rate. When the daughter branches are asymmetric in geometry but with the same mass flow rate, as shown in Figure 4, the blood cell moves to the one with a smaller cross section area, which is due to the higher average velocity in that daughter branch.

A multi-blood cell case is also simulated and the result is shown in Figure 5. The cell-to-cell interaction can be

Figure 2. Bifurcation geometry.


Figure 3. A RBC owing in a symmetric bifurcated vessel with. (a) t = 0.008 s; (b) t = 0.012 s; (c) t =0.016 s; (d) t = 0.020 s.


Figure 4. A RBC flowing in an asymmetric bifurcation vessel with. (a) t = 0.00 s; (b) t = 0.01 s; (c) t =0.02 s; (d) t = 0.03 s.


Figure 5. Two RBCs in a symmetric bifurcation. (a) t = 0.004 s; (b) t = 0.008 s; (c) t =0.012 s; (d) t = 0.018 s.

thought as two parts: 1) one cell changes the deformation and motion of the other cells by directly contacting each other; 2) one cell affects the other cells indirectly by changing the surrounding fluid field. For this case, instead of specifying the flow rate ratio between the daughter vessels, a constant inlet flow rate is given. Therefore, by simulating two RBCs laying in a line profile going through the symmetric bifurcation together, we are able to see this indirect cell-to-cell interaction. From these simulation results, the IFEM is proved to be suitable to simulate the red blood cell motion and deformation in microvessel bifurcation.

6.2. Deployment of Angioplasty Stents

For individuals with an occlusive vascular disease, blood flowing to an organ or to a distal body part is impaired by narrowed arteries with fatty deposits or calcium accumulations. Angioplasty was introduced by Dr. Andreas Gruentzig in the mid to late 1970’s and is widely used today. The area of arterial blockage is dilated with the help of a catheter that has an inflatable small balloon at its tip. Then, the plaque is squeezed along the artery wall. A decade later, stenting was introduced by Dr. Julio Palmaz in 1988 to improve the angioplasty procedure. Like angioplasty, coronary stents physically open the channel of constricted arterial segments. During stenting, a catheter delivers a balloon and a surrounding stent to the location of the blockage area. The balloon deploys the stent, remains inflated for 30 seconds and then is deflated. At the end of the process, the expanded stent is embedded into the wall of the diseased artery and holds it open.

Here, we mainly focus on the deployment process of balloon-expanding stents. Balloon expandable stents are typically made of stainless steel tubing mounted over an angioplasty balloon and then plastically deformed to their final diameter by high pressure balloon inflation.

Colombo et al. [39] and Goldberg et al. [40] demonstrated that stent apposition was inappropriate in up to 87% of the cases using conventional balloon implantation techniques. Therefore, researchers started to study the role of stent deployment in order to define an optimal deployment technique. Deviations from normal blood flow pattern may favor the initiation and progression of a vascular wall lesion. These conditions are certainly fulfilled when a stent is deployed in the vessel wall, completely or incompletely. Segers et al. [41] showed experimentally the importance of optimal stent deployment for steady state condition and pulsatile flow by studying the most important stents parameters that influence hemodynamics. Russo et al. [42] suggested that the use of high pressure non-compliant balloon stent deployment techniques give better initial results. However, these stents provoke an increased intimal growth due to more profound and deep vascular injury. Shear stress, exerted on both the vascular wall and on the blood constituents, is considered to play an important role in restenosis. Therefore, besides biocompatibility issues, the design of both the stent and delivery mechanisms should receive the necessary attention during the development process of the future generation of stents.

In our previous work [27], we built a computer model and simulated a balloon expandable stent interacting with its surrounding fluid using the immersed finite element method. The initial configuration of the balloon and stent are designed and discretized using SolidWorks [43] and imported into our IFEM code. With the previously derived computational algorithm, the displacement and stress distribution of the stent and the velocity profile of the fluid domain can be obtained and analyzed. A balloon is designed with deflated tips at its two ends in its initial undeformed configuration, as shown in Figure 6(a). A catheter is located inside the balloon to apply appropriate pressure to inflate the balloon. The balloon has a length of 15 mm, an outer diameter of 1.54 mm, and a thickness of 0.04 mm. Both ends of the balloon are fixed in all directions as shown in Figure 6(b).

The balloons used for stenting are made of very stiff polyamide (nylon) material [44,45]. The balloon is modeled as hyperelastic material with Mooney-Rivlin description in the simulation. The parameters and material properties used for the balloon are summarized in Table 1.

A stent is a cylindrical and symmetrical assembly of inter-connected diamond-shaped members. In this particular model, we use the Medtronic AVE Modular stents S7 (Medtronic AVE, Inc., Santa Rosa, CA, USA), Figure 7(a). Although this stent is no longer been widely used, it is still a good representation of a typical geometrical shape of a stent, the model can be simply modified to other expandable devices. The stent, made of wires that form a diamond shape as shown on Figure 7(b), has an outer diameter of 1.64 mm. The stent is made of 16 identical structural members with a total length of 8 mm before its expansion. The cross-section of the wire has a width of 0.08 mm. The structural members are connected peak-topeak at their tips. Balloon-expandable stents are made from materials that can be plastically deformed through the inflation of a balloon. An ideal stent should have low yield stress (to make it deformable at manageable balloon pressures), high elastic modulus (for minimal recoil), and high strength through expansion. The most widely used material for this type of stents is stainless steel, typically 316 L, a particular corrosion-resistant material. For stents the fatigue resistance is of high importance. Therefore, the microstructure requires very small

Table 1. Parameters and material properties used for the balloon.


Figure 6. Balloon geometry setup. (a) Initial balloon geometry before inflation; (b) Fixed boundaries applied at the two ends of the balloon.

(a) (b)

Figure 7. Design of a Medtronic AVE Modular stents S7.

grain sizes. The typical average grain sizes of approximately 25 μm exist in stainless steel stents and 8 to 9 grains over the strut wall thickness are desired [46]. We assume a clean homogeneous material with a uniformly fine grain size and a near 100% density. Stainless-steel alloys are usually preferred to design stents because they have been proven to be biocompatible for long-term implants in the human body and easily deformable in fully annealed condition. The stent is discretized using 9.979 nodes and 16.888 elements. A summary of the properties and parameters used for the stent is listed in Table 2.

Finally, the stent is mounted around the balloon as shown in Figure 8. The stent is placed at the center of the balloon. At this point, several assumptions must be made: 1) the force applied on the balloon is transferred directly to the stent; 2) there is no fluid between the balloon and the stent; and 3) the contact surface between the balloon and the stent is frictionless.

The fluid domain represents blood in the artery where the stent is deployed. We consider a straight artery segment with a length of 16 mm and a lumen diameter of

Table 2. Parameters and material properties used for the stent.

Figure 8. Initial configuration of catheter, ballon, and stent.

4.5 mm. Our computational model assumes that the fluid is everywhere in the domain. Although the balloon is not inflated with the same fluid as the surrounding blood, a liquid (with color) that has similar properties as the blood is normally inserted into the balloon for easier visualization in surgeries. Therefore, it is a reasonable assumption to have homogeneous fluid in the entire computational domain. Blood is an incompressible fluid consisting of a suspension of deformable particles (blood cells, platelets, etc.) in a Newtonian liquid (plasma). Except within the microcirculation, blood can be considered homogeneous with constant density and viscosity. The parameters and properties used for the fluid domain to perform the simulation can be found in Table 3.

The crimped stent is deformed by a radial force from the balloon. Experiments performed by Dumoulin and Cochelin [47], and Moore and Berry [48] showed that, except at the tips, the structure is almost uniformly dilated and finally expanded. Thus, it seems justifiable to model expansion by considering a long structure under uniform radial internal pressure. This pressure difference is applied from the centerline towards the outer diameter of the balloon and the artery wall, as shown on Figure 9. The applied inflation pressure is dependent on the struts of the wires and the specific material used. The total inflation pressure for the balloon is equal to the sum of the pressure required to deform the balloon material and to deploy the stent. During the entire simulation the inflation pressure is constant and equals to 100 g/cm2.

Our primary goal is to study the deformation of the stent during its deployment. Figure 10(a) shows the deployment of the stent during balloon expansion at different time steps. The pressure applied onto the fluid inflates the balloon and the balloon provides a force onto the stent, which enables the stent to expand radially outward until it contacts the inner surface of the artery wall.

During deployment, the diameter of the stent increases from 1.64 mm to 2.82 mm. As expected, the stent deforms uniformly in its radial direction, expanding by 1.7 times its initial diameter for our given material. Figure 11(a) shows that the expansion of the stent in its radial direction follows a linear variation.

Figures 10(a) demonstrates how the diamond shape of

Table 3. Parameters and properties used for the fluid domain.

Figure 9. The initial conditions applied to the fluid domain: constant pressure difference (P2 − P1) is applied from the catheter to the fluid boundaries.


Figure 10. Stent configuration and von Mises stress map during stent deployment. (a) Stent deployment process; (b) Von Mises stress during stent deployment.


Figure 11. Stent diameter and maximum von Mises stress vs. time during stent deployment. (a) Stent diameter variation during the deployment; (b) Maximum Von Mises stress variation during the deployment.

the stent deforms during the deployment. Compared to experiments and observations of deployment of a real stent, our three-dimensional computational model gives a realistic representation of stent expansion during balloon inflation.

The strength and the long term in-vivo performance of the stent can be determined from the stress distributions with the goal of minimizing vascular injury. The Von Mises stress distribution along the stent is reported at different time steps in Figure 10(b). The stress distribution is uniform longitudinally along the stent and varies during deployment. The highest stress values appear at the final stage of expansion, and the peak Von Misses stress follows a linear variation during the entire simulation as shown on these values are critical for recoil and failure analysis. It can be concluded that the stress distribution is a function of applied pressure, balloon and stent material properties, fluid properties, and stent geometry. Figure 12 shows the radial fluid velocity profile during the deployment of the stent. This profile represents better how the applied uniform pressure acts on the balloon in the radial direction. Figure 13 shows the magnitude and profile of the velocity in the transverse plane. Both figures demonstrate a uniform velocity distribution in the direction of expansion, illustrating how well our model simulates the expansion mechanism through applied pressure difference.

Figure 12. Radial fluid velocity profile at different time steps.

Figure 13. Velocity profile in the longitudinal direction during expansion at different time steps.

6.3. Human Vocal Folds Vibration during Phonation

Voice is produced by the vibration of vocal folds. The vocal folds are a pair of pliable structures located within the larynx at the top of the trachea, see Figure 14. The

Figure 14. Human vocal folds.

human vocal folds are roughly 10 - 15 mm in length and 3 - 5 mm thick. The human vocal folds are laminated structures composed of five different layers: the epithelium, the superficial layer (SLP), the intermediate layer (ILP), the deep layer and thyroarytenoid muscle, shown in Figure 14.

An accurate numerical simulation of the vocal folds vibration can help us obtain a better understanding of the dynamics of the voice production in human beings. Due to the complicated nature of this problem, the numerical model has to fulfill the following requirements. First, the numerical model has to represent completely a coupled fluid-structure interaction system. Second, the numerical model should perform well when there exists large density ratio between the fluid and structure because the density of the vocal fold muscle is close to water and the density ratio between the vocal fold muscle and the airflow is about 1000. Third, the motion and deformation of the structure have to be predicted accurately with complicated geometry and material descriptions because the vocal folds have complex shape and layer-structures and are viscoelastic materials. The mIFEM is a perfect numerical method to perform the simulation of this complex problem.

The geometry of the self-oscillated vocal folds model is shown in Figure 15. Since sound is generated by the compression of air, the working fluid is taken as compressible air governed by the ideal gas law at room temperature. The density of the fluid is and the viscosity of the fluid is.

The vocal fold muscle is considered as isotropic viscoelastic material. The vocal fold is assumed to have layered structure, outside cover layer (red) and inside body layer (green). The cover layer is much softer than the body layer. For the cover layer the Young’s modulus is, whereas for the body layer. The densities of both cover and body layer are assumed to be the same as. The Poisson ratio is. Two vocal folds have the exact same geometry

Figure 15. 2-D two-layer self-oscillated vocal folds model.

and material description, sit in the fluid channel symmetric about the central line. A constant total pressure boundary condition of is applied at the channel inlet and the outflow boundary is given at the channel exit. No-slip and no-penetration boundary conditions are applied on the channel walls and on the vocal fold surfaces.

A snapshot of the fluid velocity field at two typical instances during a steady vibration cycle are shown in Figure 16. One can see that the fluid field is not symmetric about the central line during the vibration. The glottal jet tends to attach to one side of the vocal folds randomly, which is the so-called the “Coanda effect” [49,50].

The asymmetrical airflow causes an asymmetrical pressure distribution in regions near the vocal folds and change in the vibration pattern. The minimum distance between the vocal fold surface and the central line is measured to represent the half glottis width, shown in Figure 17, where Gwup and Gwdown represent the opening width for the up and down vocal folds, respectively. This figure shows that the simulation captures the vocal folds to have a repeated opening and closing process. When the glottis width is zero or near-zero, then the vocal folds are closed, there is no air flowing through. The pressure starts to build up in the upstream of the vocal channel. As the pressure increases, it starts to push the vocal folds to open and eventually reaches a maximum glottis width, the high pressure is released. The vocal folds then return back to its closed position, and the whole process restarts. The glottis width for the up and down vocal folds do not equal each other over cycles, indicating that the vocal folds motion is asymmetric although the vibrational frequency is found to be the same. The vibrational magnitude is slightly off from each other. To find out the vibration frequency, FFT is performed on the up and down glottis widths and the volume flow rate,. The power spectra are plotted in Figure 18 where the frequency for these three variables are the same and found to be 234 Hz, which is in the expected range of a female vocal fold vibrational frequency.


In this paper, we reviewed IFEM algorithms that had

Figure 16. Fluid velocity field at two typical instances during steady vibration.

Figure 17. Half glottis width of top and bottom vocal folds.

Figure 18. Spectrum plot of flow rate and half glottis width of the top and bottom vocal folds.

been developed over the past decade. The IFEM method is a numerical scheme that adopts the non-boundary-fitted mesh approach and fully couples the fluid-structure interaction by interpolation of the interacting domains. The fluid and solid domains are solved independently using finite element method and coupled with each other within one time step through fluid-structure interaction force. The original IFEM algorithm is considered as explicitly coupling for the fluid and solid, which can lead to instability issue when time step is not sufficiently small. The semi-implicit IFEM algorithm tackles this problem by modifying the FSI force and adopting the concept of the indicator function. It extends the stability range of the numerical scheme and allows us to consider the fluidstructure interaction problems when the fluid and solid properties are very different from each other, for example high density ratio between the solid and fluid, and the solid with relatively large stiffness. For high Reynolds number flows,, the modified IFEM algorithm performs better due to its accurate description of the solid motion and deformation through capturing the dynamics of the solid motion. Three biomedical applications are studied using these IFEM algorithms to simulate soft tissues interacting with fluid where the fluid can be either blood or air.


This work is partially funded by NIH (R01-DC005642-07).


  1. Peskin, C. (1977) Numerical analysis of blood flow in the heart. Journal of Computational Physics, 25, 220-252.
  2. McCracken, M. and Peskin, C. (1980) A vortex method for blood flow through heart valves. Journal of Computational Physics, 35, 183-205.
  3. McQueen, D. and Peskin, C. (1983) Computer-assisted design of pivoting-disc prosthetic mitral valves. Journal of Computational Physics, 86, 126-135.
  4. Peskin, C. and McQueen, D. (1989) A three-dimensional computational method for blood flow in the heart. I. Immersed elastic fibers in a viscous incompressible fluid. Journal of Computational Physics, 81, 372-405.
  5. Peskin, C. and McQueen, D. (1992) Cardiac fluid dynamics. Critical Reviews in Biomedical Engineering, SIAM Journal on Scientific and Statistical Computing, 20, 451-459.
  6. Peskin, C. and McQueen, D. (1994) Mechanical equilibrium determines the fractal fiber architecture of aortic heart valve leaflets. American Journal of Physiology, 266, H319-H328.
  7. Peskin, C. and McQueen, D. (1996) Case studies in mathematical modeling-ecology, physiology, and cell biology. Prentice-Hall, Upper Saddle River.
  8. LeVeque, R.J. and Li, Z. (1994) The immersed interface method for elliptic equations with discontinuous coefficients and singular sources. SIAM Journal on Numerical Analysis, 31, 1091-1044.
  9. LeVeque, R.J. and Li, Z. (1997) Immersed interface methods for stokes flow with elastic boundries or surface tension. SIAM Journal on Scientific Computing, 18, 709- 735.
  10. Fogelson, A. and Keener, J. (2000) Immersed interface method for Neumann and related problems in two and three dimensions. SIAM Journal on Scientific Computing, 22, 1630-1654.
  11. Lee, L. and LeVeque, R. (2003) An immersed interface method for incompressible Navier-Stokes equations, SIAM Journal on Scientific Computing, 25, 832-856.
  12. Li, Z. and Lai, M. (2001) The immersed interface method for the Navier-Stokes equations with singular forces. Journal of Computational Physics, 171, 822-842.
  13. Wiegmann, A. and Bube, K.P. (1998) The immersed interface method for nonlinear differential equations with discontinuous coefficients and singular sources. SIAM Journal on Numerical Analysis, 35, 177-200.
  14. Wiegmann, A. and Bube, K.P. (2000) The explicit-jump immersed interface method: Finite difference methods for PDEs with piecewise smooth solutions. SIAM Journal on Numerical Analysis, 37, 827-862.
  15. Wang, X. and Liu, W. (2004) Extended immersed boundary method using FEM and RKPM. Computer Methods in Applied Mechanics and Engineering, 193, 1305-1321.
  16. Boffi, D. and Gastaldi, L. (2003) A finite element approach for the immersed boundary method. Computers and Structures, 81, 491-501.
  17. Boffi, D., Gastaldi, L. and Heltai, L. (2007) On the CFL condition for the finite element immersed boundary method. Computers and Structures, 85, 775-783.
  18. Wang, S.X., Zhang, L.T. and Liu, W.K. (2009) Finite element formulations for immersed methods: Explicit and implicit approaches. Journal of Computational Physics, 228, 2535-2551.
  19. Zhang, L., Gerstenberger, A., Wang, X. and Liu, W. (2004) Immersed finite element method. Computer Methods in Applied Mechanics and Engineering, 193, 2051-2067.
  20. Zhang, L. and Gay, M. (2007) Immersed finite element method for fluid-structure interactions. Journal of Fluids and Structures, 23, 839-857.
  21. Zhang, L.T. and Gay, M. (2008) Imposing rigidity constraints on immersed objects in unsteady fluid flows. Computational Mechanics, 42, 357-370.
  22. Wang, X. and Zhang, L.T. (2010) Interpolation functions in the immersed boundary and finite element methods. Computational Mechanics, 45, 321-334.
  23. Liu, W., Liu, Y., Farrell, D., Zhang, L., Wang, S., Fukui, Y., Patankar, N., Zhang, Y., Bajaj, C., Lee, J., Hong, J., Chen, X. and Hsu, H. (2006) Immersed finite element method and its applications to biological systems. Computer Methods in Applied Mechanics and Engineering, 195, 1722-1749.
  24. Liu, W., Liu, Y., Zhang, L., Wang, X., Gerstenberger, A. and Farrell, D. (2004) Immersed finite element method and applications to biological systems. Finite Element Methods: 1970’s and Beyond. International Center for Numerical Methods and Engineering.
  25. Liu, Y. and Liu, W. (2006) Rheology of red blood cell aggregation in capillary by computer simulation. Journal of Computational Physics, 220, 139-154.
  26. Liu, Y., Zhang, L., Wang, X. and Liu, W. (2004) Coupling of Navier-Stokes equations with protein molecular dynamics and its application to hemodynamics. International Journal for Numerical Methods in Fluids, 46, 1237-1252.
  27. Gay, M., Zhang, L. and Liu, W. (2006) Stent modeling using immersed finite element method. Computer Methods in Applied Mechanics and Engineering, 195, 4358- 4370.
  28. Zhang, L.T. (2008) Shear stress and shear-induced particle residence in stenosed blood vessels, International Journal of Multiscale Computational Engineering, 6, 141-152.
  29. Zhang, L.T. and Gay, M. (2008) Characterizing left atrial appendage functions in sinus rhythm and atrial fibrillation using computational models. Journal of Biomechanics, 41, 2515-2523.
  30. Gay, M. and Zhang, L.T. (2009) Numerical studies of healthy, stenosed, and stented coronary arteries. International Journal of Numerical Methods in Fluids, 61, 453- 472.
  31. M. Gay, and L. T. Zhang, (2009) Numerical studies on fluid-structure interactions of stent deployment and stented arteries. Engineering with Computers, 25, 61-72.
  32. Tezduyar, T. (1992) Stabilized finite element formulations for incompressible-flow computations. Advanced Application Mechanics, 28, 1-44.
  33. Tezduyar, T. (2001) Finite element methods for flow problems with moving boundaries and interfaces. Archives of Computational Methods in Engineering, 8, 83-130.
  34. Hughes, T., Franca, L. and Balestra, M. (1986) A new finite element formulation for computational fluid dynamics: V. Circumventing the babuška-brezzi condition: A stable Petrov-Galerkin formulation of the stokes problem accommodating equal-order interpolations. Computer Methods in Applied Mechanics and Engineering, 59, 85-99.
  35. Wang, X., Wang, C. and Zhang, L.T. (2011) Semi-implicit formulation of the immersed finite element method. Computational Mechanics, 49, 421-430.
  36. Wang, X. and Zhang, L.T. (2013) Modified immersed finite element method for solid-dominated fully-coupled fluid-structure interactions. Computer Methods in Computer Methods in Applied Mechanics and Engineering, 267, 150-169.
  37. Peskin, C. (2002) The immersed boundary method. Acta Numerica, 11, 479-517.
  38. Torres, D. and Brackbill, J. (2000) The point-set method: front-tracking without connectivity. Journal of Computational Physics, 165, 620-644.
  39. Colombo, A., Hall, P., Nakamura, S., Almagor, Y., Maiello, L., Martini, G., Gaglione, A., Goldberg, S. and Tobis, J. (1995) Intracoronary stenting without anticoagulation accomplished with intravascular ultrasound guidance. Circulation, 91, 1676-1688.
  40. Goldberg, S., Colombo, A., Nakamura, S., Almagor, Y., Maiello, L. and Tobis, J. (1994) Benefit of intracoronary ultrasound in the deployment of Palmaz-Schatz stents. Journal of the American College of Cardiology, 24, 996- 1003.
  41. Segers, P., Kostopoulos, K., Scheerder, I.D. and Verdonck, P. (1998) Biomechanical aspects of intracoronary stents. In: Verdonck, P., Ed., Intra and Extracorporeal Cardiovascular Fluid Dynamics, Vol. 1, General Principles in Application, WIT Press, Ashurst, 203-232.
  42. Russo, R., Schatz, R., Sklar, M., Johnson, A., Tobis, J. and Teirstein, P. (1995) Ultrasound guided coronary stent placement without prolonged systemic anticoagulation. Journal of the American College of Cardiology, 25, 50A.
  43. SolidWorks (2004) Sp03.1. SolidWorks Corporation, Concord.
  44. Serruys, P. and Rensing, B. (2002) Handbook of coronary stents. 4th Edition, Martin Dunitz, London.
  45. Saab, M. (1999) Applications of high-pressure balloons in the medical device industry. Advanced Polymers, Inc. Salem.
  46. Murphy, B., Savage, P., McHugh, P. and Quinn, D. (2003) The stress-strain behavior of coronary stent struts is size dependent. Annals of Biomedical Engineering, 31, 686- 691.
  47. Dumoulin, C. and Cochelin, B. (2000) Mechanical behavior modeling of balloon-expandable stents. Journal of Biomechanics, 33, 1461-1470.
  48. Moore Jr., J.E. and Berry, J. (2002) Fluid and solid mechanical implications of vascular stenting. Annals of Biomedical Engineering, 30, 498-508.
  49. Tao, C., Zhang, Y., Hottinger, D.G. and Jiang, J.J. (2007) Asymmetric airflow and vibration induced by the Coanda effect in a symmetric model of the vocal folds. Journal of the Acoustical Society of America, 122, 2270-2278.
  50. Drechsel, J.S. and Thomson, S.L. (2008) Influence of supraglottal structures on the glottal jet exiting a twolayer synthetic, self-oscillating vocal fold model. Journal of the Acoustical Society of America, 123, 4434-4445.