The interactions between incompressible fluid flows and immersed structures are nonlinear multi-physics phenomena that have applications to a wide range of scientific and engineering disciplines. In this article, we review representative numerical methods based on conforming and non-conforming meshes that are currently available for computing fluid-structure interaction problems, with an emphasis on some of the recent developments in the field. A goal is to categorize the selected methods and assess their accuracy and efficiency. We discuss challenges faced by researchers in this field, and we emphasize the importance of interdisciplinary effort for advancing the study in fluid-structure interactions.

Explicit time stepping schemes for the immersed boundary method require very small time steps in order to maintain stability. Solving the equations that arise from an implicit discretization is difficult. Recently, several different approaches have been proposed, but a complete understanding of this problem is still emerging. A multigrid method is developed and explored for solving the equations in an implicit-time discretization of a model of the immersed boundary equations. The model problem consists of a scalar Poisson equation with conformation-dependent singular forces on an immersed boundary. This model does not include the inertial terms or the incompressibility constraint. The method is more efficient than an explicit method, but the efficiency gain is limited. The multigrid method alone may not be an effective solver, but when used as a preconditioner for Krylov methods, the speed-up over the explicit-time method is substantial. For example, depending on the constitutive law for the boundary force, with a time step 100 times larger than the explicit method, the implicit method is about 15-100 times more efficient than the explicit method. A very attractive feature of this method is that the efficiency of the multigrid preconditioned Krylov solver is shown to be independent of the number of immersed boundary points.

The immersed boundary (IB) method is an approach to problems of fluid-structure interaction in which an elastic structure is immersed in a viscous incompressible fluid. The IB formulation of such problems uses a Lagrangian description of the structure and an Eulerian description of the fluid. It is well known that some versions of the IB method can suffer from poor volume conservation. Methods have been introduced to improve the volume-conservation properties of the IB method, but they either have been fairly specialized, or have used complex, nonstandard Eulerian finite-difference discretizations. In this paper, we use quasi-static and dynamic benchmark problems to investigate the effect of the choice of Eulerian discretization on the volume-conservation properties of a formally second-order accurate IB method. We consider both collocated and staggered-grid discretization methods. For the tests considered herein, the staggered-grid IB scheme generally yields at least a modest improvement in volume conservation when compared to cell-centered methods, and in many cases considered in this work, the spurious volume changes exhibited by the staggered-grid IB method are more than an order of magnitude smaller than those of the collocated schemes. We also compare the performance of cell-centered schemes that use either exact or approximate projection methods. We find that the volume-conservation properties of approximate projection IB methods depend strongly on the formulation of the projection method. When used with the IB method, we find that pressure-free approximate projection methods can yield extremely poor volume conservation, whereas pressure-increment approximate projection methods yield volume conservation that is nearly identical to that of a cell-centered exact projection method.

Many problems involving the interaction of an elastic structure and a viscous fluid can be solved by the immersed boundary (IB) method. In the IB approach to such problems, the elastic forces generated by the immersed structure are applied to the surrounding fluid, and the motion of the immersed structure is determined by the local motion of the fluid. Recently, the IB method has been extended to treat more general elasticity models that include both positional and rotational degrees of freedom. For such models, force and torque must both be applied to the fluid. The positional degrees of freedom of the immersed structure move according to the local linear velocity of the fluid, whereas the rotational degrees of freedom move according to the local angular velocity. This paper introduces a spatially adaptive, formally second-order accurate version of this generalized immersed boundary method. We use this adaptive scheme to simulate the dynamics of an elastic ring immersed in fluid. To describe the elasticity of the ring, we use an unconstrained version of Kirchhoff rod theory. We demonstrate empirically that our numerical scheme yields essentially second-order convergence rates when applied to such problems. We also study dynamical instabilities of such fluid-structure systems, and we compare numerical results produced by our method to classical analytic results from elastic rod theory.

The immersed boundary method has been extensively used to simulate the motion of elastic structures immersed in a viscous fluid. For some applications, such as modeling biological materials, capturing internal boundary viscosity is important. We present numerical methods for simulating Kelvin-Voigt and standard linear viscoelastic structures immersed in zero Reynolds number flow. We find that the explicit time immersed boundary update is unconditionally unstable above a critical boundary to fluid viscosity ratio for a Kelvin-Voigt material. We also show there is a severe time step restriction when simulating a standard linear boundary with a small relaxation time scale using the same explicit update. A stable implicit method is presented to overcome these computation challenges.

We extend the immersed boundary (IB) method to simulate the dynamics of a 2D dry foam by including the topological changes of the bubble network. In the article [Y. Kim, M.-C. Lai, and C. S. Peskin, J. Comput. Phys. 229: 5194-5207, 2010], we implemented an IB method for the foam problem in the two-dimensional case, and tested it by verifying the von Neumann relation which governs the coarsening of a two-dimensional dry foam. However, the method implemented in that article had an important limitation; we did not allow for the resolution of quadruple or higher order junctions into triple junctions. A total shrinkage of a bubble with more than four edges generates a quadruple or higher order junction. In reality, a higher order junction is unstable and resolves itself into triple junctions. We here extend the methodology previously introduced by allowing topological changes, and we illustrate the significance of such topological changes by comparing the behaviors of foams in which topological changes are allowed to those in which they are not.

We present a three dimensional model of an open elastic tube immersed in fluid to understand valveless pumping mechanism. A fluid-tube interaction problem is simulated by the volume conserved immersed boundary method which prevents the generation of spurious velocity field near the tube and local cluster of the tube surface. In order to explain pumping phenomena without valves, average net flow is measured by changing parameter values such as pumping frequency, compression duration, and pumping amplitude. Some frequencies that make the system reach maximal or minimal net flow are selected to study case by case. We also study the effectiveness of fluid mixing using the Shannon entropy increase rate.

An adaptive mesh refinement strategy is proposed in this paper for the Immersed Boundary and Immersed Interface methods for two-dimensional elliptic interface problems involving singular sources. The interface is represented by the zero level set of a Lipschitz function $ϕ(x,y)$. Our adaptive mesh refinement is done within a small tube of $|ϕ(x,y)|≤δ$ with finer Cartesian meshes. The discrete linear system of equations is solved by a multigrid solver. The AMR methods could obtain solutions with accuracy that is similar to those on a uniform fine grid by distributing the mesh more economically, therefore, reduce the size of the linear system of the equations. Numerical examples presented show the efficiency of the grid refinement strategy.

An iterative solver based on the immersed interface method is proposed to solve the pressure in a two-fluid flow on a Cartesian grid with second-order accuracy in the infinity norm. The iteration is constructed by introducing an unsteady term in the pressure Poisson equation. In each iteration step, a Helmholtz equation is solved on the Cartesian grid using FFT. The combination of the iteration and the immersed interface method enables the solver to handle various jump conditions across two-fluid interfaces. This solver can also be used to solve Poisson equations on irregular domains. 

A novel full Eulerian fluid-elastic membrane coupling method on the fixed Cartesian coordinate mesh is proposed within the framework of the volume-of-fluid approach. The present method is based on a full Eulerian fluid-(bulk) structure coupling solver (Sugiyama et al., J. Comput. Phys., 230 (2011) 596–627), with the bulk structure replaced by elastic membranes. In this study, a closed membrane is considered, and it is described by a volume-of-fluid or volume-fraction information generally called VOF function. A smoothed indicator (or characteristic) function is introduced as a phase indicator which results in a smoothed VOF function. This smoothed VOF function uses a smoothed delta function, and it enables a membrane singular force to be incorporated into a mixture momentum equation. In order to deal with a membrane deformation on the Eulerian mesh, a deformation tensor is introduced and updated within a compactly supported region near the interface. Both the neo-Hookean and the Skalak models are employed in the numerical simulations. A smoothed (and less dissipative) interface capturing method is employed for the advection of the VOF function and the quantities defined on the membrane. The stability restriction due to membrane stiffness is relaxed by using a quasi-implicit approach. The present method is validated by using the spherical membrane deformation problems, and is applied to a pressure-driven flow with the biconcave membrane capsules (red blood cells).

In this paper, we investigate the dynamic process of liquid bridge formation between two parallel hydrophobic plates with hydrophilic patches, previously studied in [1]. We propose a dynamic Hele-Shaw model to take advantage of the small aspect ratio between the gap width and the plate size. A constrained level set method is applied to solve the model equations numerically, where a global constraint is imposed in the evolution [2] stage together with local constraints in the reinitialization [3] stage of level set function in order to limit numerical mass loss. In contrast to the finite element method used in [2], we use a finite difference method with a 5th order HJWENO scheme for spatial discretization. To illustrate the effectiveness of the constrained method, we have compared the results obtained by the standard level set method with those from the constrained version. Our results show that the constrained level set method produces physically reasonable results while that of the standard method is less reliable. Our numerical results also show that the dynamic nature of the flow plays an important role in the process of liquid bridge formation and criteria based on static energy minimization approach has limited applicability.

Solving elasticity equations with interfaces is a challenging problem for most existing methods. Nonetheless, it has wide applications in engineering and science. An accurate and efficient method is desired. In this paper, an efficient non-traditional finite element method with non-body-fitting grids is proposed to solve elasticity equations with interfaces. The main idea is to choose the test function basis to be the standard finite element basis independent of the interface and to choose the solution basis to be piecewise linear satisfying the jump conditions across the interface. The resulting linear system of equations is shown to be positive definite under certain assumptions. Numerical experiments show that this method is second order accurate in the L^∞ norm for piecewise smooth solutions. More than 1.5th order accuracy is observed for solution with singularity (second derivative blows up) on the sharp-edged interface corner.