A stable fluid-structure-interaction solver for low-density rigid bodies using the immersed boundary projection method

# A stable fluid-structure-interaction solver for low-density rigid bodies using the immersed boundary projection method

[    [    [
###### Abstract

Dispersion of low-density rigid particles with complex geometries is ubiquitous in both natural and industrial environments. We show that while explicit methods for coupling the incompressible Navier-Stokes equations and Newton’s equations of motion are often sufficient to solve for the motion of cylindrical particles with low density ratios, for more complex particles – such as a body with a protrusion – they become unstable. We present an implicit formulation of the coupling between rigid body dynamics and fluid dynamics within the framework of the immersed boundary projection method. Similarly to previous work on this method, the resulting matrix equation in the present approach is solved using a block-LU decomposition. Each step of the block-LU decomposition is modified to incorporate the rigid body dynamics. We show that our method achieves second-order accuracy in space and first-order in time (third-order for practical settings), only with a small additional computational cost to the original method. Our implicit coupling yields stable solution for density ratios as low as . We also consider the influence of fictitious fluid located inside the rigid bodies on the accuracy and stability of our method.

$\ast$$\ast$footnotetext: Corresponding author. Tel.: +46 8 790 71 77

A stable fluid-structure-interaction solver for low-density rigid bodies using the immersed boundary projection method

aLinné Flow Centre, KTH Mechanics, 10044 Stockholm, Sweden

bDepartment of Mechanical Engineering and Florida Center for Advanced Aero-Propulsion, Florida State University, Tallahassee, Florida 32310, USA

Keywords: Immersed boundary method, Fictitious fluid, Newton’s equations of motion, Implicit coupling, Low density ratios, Complex particles

## 1 Introduction

During recent years, the original immersed boundary (IB) method peskin1972flow () has been extensively developed and gained popularity due to the ability to handle the interaction of objects of complex geometries with fluids. The key feature of the IB method is that the underlying Eulerian grid does not need to be body conforming. The application of the method ranges from fundamental problems of solid particle suspensions Uhlmann_JCP_2005 (); kempe2012improved (); Breugem_JCP_2012 (), to natural and industrial problems of complex and elastic geometries peskin1977numerical (); ye1999accurate (); kim2001immersed (); tseng2003ghost (); griffith2012immersed (); borazjani2013fluid ().

A number of studies Causin20054506 (); Forster20071278 (); Borazjani20087587 () have reported difficulties with numerical convergence in fluid-structure interaction (FSI) problems, when the fluid force acting on the solid dominates over the solid inertial force, for example, in blood flow through flexible arteries, as discussed by Baek and Karniadakis Baek2012629 (). Often the numerical instabilities are attributed to the added mass component. In the IB framework convergence problems have been reported by Borazjani et al. Borazjani20087587 ().

Numerical instabilities are also present in particulate flows. Uhlman Uhlmann_JCP_2005 () developed an efficient direct forcing IB method that describes the coupling between the Newton’s equations of motion and the incompressible Navier-Stokes equations for many particles. However, due to the assumption that the fluid inside the particle moves as a solid body, the method suffers from stability issues for density ratio between solid and fluid below (for spheres), where is the density of the solid and is the density of the fluid. The method was later improved by Kempe and Fröhlich kempe2012improved () by taking the motion of the fluid inside the particle explicitly into account; they were able to reduce the critical density ratio to . Same approach was also used by Breugem Breugem_JCP_2012 (). All of these approaches use explicit coupling (so called weak coupling) between the rigid body and the fluid. The works by Uhlman Uhlmann_JCP_2005 (), Kempe and Fröhlich kempe2012improved () and Breugem Breugem_JCP_2012 () do not consider non-spherical bodies for which the fluid forces acting on the solid can be significantly larger. For such bodies, their algorithms are likely to be stable only for heavier particles with higher .

Some studies have investigated the dynamics of more complex (non-spherical) particle geometries Borazjani20087587 (); zheng2010coupled (); yang2015non (); Gibou20123246 (). Zheng et al. zheng2010coupled () investigated human phonation numerically. To solve the motion of human vocal cords and interaction with surrounding air, they have developed a numerical method, in which a sharp-interface IB method is explicitly coupled with-finite difference Navier-Stokes solver. They show that the limiting density ratio, for which the method becomes unstable, is , which is similar as in the work by Kempe and Fröhlich kempe2012improved (). They also have derived the necessary time step in order to resolve the motion of vocal cords and preserve stability. Borazjani et al. Borazjani20087587 () compared weak coupling (WC) and strong coupling (SC) algorithms within the IB framework for different problems. The SC in their method is ensured using Gauss-Seidel-like iterations within each time step. They noticed that the WC algorithm becomes unstable, when the mass of the solid structure is reduced below some critical value. For certain problems they noticed that the SC algorithm also suffers from a similar drawback. A relaxation scheme for the inner iterations was implemented to overcome the problem, but with increased number of iterations and added computational cost.

Yang and Stern yang2015non () have very recently presented a strongly-coupled non-iterative method using fractional step approach to solve Navier-Stokes equations coupled with sharp-interface direct-forcing IB method. They implemented a SC algorithm by introducing an intermediate step in a non-inertial reference frame, following the motion of solid body. Yang and Stern demonstrate improved stability properties with stable simulations down to density ratio , which is a significant improvement over work by Kempe and Fröhlich kempe2012improved (). The derived formulation does not require any iterations within each time step and thus reduces the computational cost. Another non-iterative method is proposed by Gibou and Min Gibou20123246 (), which is similar to algorithm described in the current paper. They advance both fluid and solid through intermediate states, and impose the interaction between fluid and solid during the projection step. In their work, the solution steps of the projection method for rigid body are deduced empirically by mimicking the fluid solution steps, while we use a more rigorous approach to derive solution steps for rigid body dynamics by a block-LU decomposition.

An alternative way to achieve better stability properties in simulations at low particle densities is including some information of added mass in the computational method. For example, Eldredge eldredge2008dynamically () has shown that the FSI computation can be stabilized if the added mass matrix is computed explicitly and added to the body inertia. Furthermore, Wang and Eldredge wang2015strongly () have used some information about the added mass to arrive with relaxation factor, which leads to stable simulations.

In the present method, we discretize the system of equations following the approach by Taira and Colonius Taira_JCP_2007 () and form a discrete linear system of equations. We then decompose the system using a block-LU decomposition, which gives us the prediction step for both fluid and solid body motion, the modified Poisson equation for the dynamic interaction force between fluid and solid, and the projection step for enforcing the interaction of the solid and fluid.

Wang and Eldredge wang2015strongly () have developed a numerical method in which the null-space fluid solver of Colonius and Taira (colonius2008fast, ) is iteratively coupled with general equations for rigid body dynamics. The null-space based IB method (colonius2008fast, ) is an extension to the original IB projection method Taira_JCP_2007 (). Our present FSI solver employs s direct solver for a positive-definite algebraic system, based on the block-LU decomposition in line with the original fractional step method Perot_JCP_1993 (), thus eliminating the need for any iterations within a single time step. This approach is illustrated using a special case of rigid body dynamics (non-deformable objects), while allowing extension to deformable, infinitely thin, open filaments and sheets. In addition, the current method gives direct access to the pressure field, which is useful in many applications.

We characterize the stability properties of our method on a vortex-induced-vibration (VIV) problem for two particles – a circular cylinder with and without a splitter plate clamped to the rear end. The current method is shown to be stable for solving the flow and body dynamics for both bodies for density ratios as low as .

In section 2, we discuss the general governing equations for the physical problem of interest. In section 3, we describe the basic elements of the IB projection method by Taira and Colonius Taira_JCP_2007 (). In section 4, we present our extension to the IB projection method for FSI problems with rigid bodies. We discuss Newton’s equations of motion and couple them with the IB projection method using both explicit (WC) and implicit formulation (SC). We formulate the SC scheme in matrix form and decompose it using a block-LU decomposition. In section 5, we show the convergence properties of the current method. We also present results of a freely falling and rising circular cylinder and a neutrally buoyant circular cylinder in shear flow and compare our findings with literature. In section 6, we demonstrate the stability properties of our method for the VIV problem mentioned previously. In section 7, we investigate the effect of fictitious fluid inside the particle on numerical stability. Finally, we draw conclusions in section 8. In “Appendix. Stability of the implicit coupling for massless particles”, we present a modification of the present method, which is stable for limiting case of massless particles. In “Appendix. Designing a parallel Poisson solver by using the block-LU decomposition”, we show a design of a parallel algorithm for solution of Poisson equation, which does not depend on domain decomposition.

## 2 Governing equations of a rigid body motion in a fluid

We consider a general solid body (represented with gray body in Fig. 1) immersed in a viscous, incompressible fluid. We denote the fluid domain with , the outer boundary with and the boundary at solid body with . The whole system can be subject to gravitational acceleration111The formulation holds also for uniform background acceleration. . The solid body moves with velocity under the influence of gravity and contact forces from the fluid.

This configuration is governed by a system of non-dimensional Navier-Stokes equations and Newton’s equations of motion

 ∂u∂t+(u⋅∇)u =−∇~p+1{Re}∇2u+g^eg in Ω, (1) ∇⋅u =0 in Ω, (2) u =us+ωs×r on S, (3) au+b∂u∂^n =c on ∂Ω, (4) dusdt =1ρVs∮S¯¯¯¯¯¯¯τ⋅^ndS+(1−1ρ)g^eg, (5) d(¯¯¯¯Is⋅ωs)dt =1ρ∮Sr×(¯¯¯¯¯¯¯τ⋅^n)dS, (6)

where equations (34) are the boundary conditions on the physical and computational domains, and , and are known parameters. The vector field is the fluid velocity field, the scalar field is the fluid pressure, (defined later) is the Reynolds number, and is the fluid stress tensor. In Newton’s equations of motion, is the translation velocity located at the center of mass for the given body, is the angular velocity of the center of mass, is the radius from the center of mass to the surface of the body or some other point in fluid, is dimensionless volume, and is dimensionless second-rank moment of inertia tensor with components .

## 3 Immersed boundary projection method

We start with writing the incompressible Navier-Stokes equations (1-2) with boundary conditions (34) in IB formulation. The solid body is replaced with fluid and volume forcing, to mimic the boundary conditions from the solid body. The system of equations is given by

 ∂u∂t+(u⋅∇)u =−∇p+1{Re}∇2u+∫SF(L)δ(L−x)dS in Ω, (7) ∇⋅u =0 in Ω, (8) u(L) =∫Ωu(x)δ(x−L)dV=us+ωs×~L, (9) u =0 on ∂Ω, (10)

where is the surface Lagrangian force density of the IB method (applied only on solid body surface or Lagrangian points), is a vector directed to the IB Lagrangian points, is a vector directed to the IB Lagrangian points from the center of solid body, is the coordinate vector in Eulerian space, and is the Dirac delta function. The effect of gravity in the Navier-Stokes equations from now on is incorporated in the modified pressure field . Without the loss of generality, we apply no-slip boundary condition (, ) on . The fluid solver used in the current work is the finite volume IB projection method on a non-uniform staggered grid ferziger2002computational (). The delta function is approximated using -cell discrete delta function developed by Roma et al. roma1999adaptive (). The diffusive term is integrated with the second-order implicit Crank–Nicolson method; the non-linear advective term is advanced in time with the second-order explicit Adams–Bashforth method. Consequently, equations (710) become

 un+1−unΔt+[32^N(un)−12^N(un−1)] =−^Gϕn+1/2+12{Re}^L(un+1+un)+ +^Hnfn+1/2+^bc1, (11) ^Dun+1 =0+^bc2, (12) ^Enun+1 =uns+ωns×~Ln, (13)

where is the non-linear (or advection) operator, is the gradient operator, is the Laplace operator, is the divergence operator, is the spreading (regularization) operator, is the interpolation operator, is the discrete pressure, is the discrete IB forcing vector, and and are boundary conditions associated with the momentum and continuity equations, respectively. Note that we have introduced the subscript to denote the time level used by the interpolation and spreading operators. This illustrates that in the case of a moving body we do not know the position of Lagrangian points a priori, therefore one option is to use positions at previous time step . In such way, we obtain a time-lagged interpolation, where the interpolation operator at time level acts on the velocity flux at time level . We will discuss this issue in further details in section 4.

The divergence and gradient matrices can be made to consist of only and by introducing diagonal scaling matrices and , i.e. and (), as detailed in Taira and Colonius Taira_JCP_2007 (). Using a similar transformation, the interpolation matrix and spreading matrix can be made transpose of each other. We introduce an implicit operator , which contains the terms in front of the unknown velocity field at the time level . Adopting the scaling above, equations (1113) of the IB formulation can be written in algebraic form as

 ⎛⎜⎝AGETnGT00En00⎞⎟⎠⎛⎜ ⎜⎝qn+1ϕn+1/2~fn+1/2⎞⎟ ⎟⎠=⎛⎜⎝rn0uns⎞⎟⎠+⎛⎜⎝bc1−bc20⎞⎟⎠, (14)

where is the velocity flux vector, is a rescaled IB forcing vector and is prescribed velocity of IB Lagrangian points. The IB forcing has been rescaled to obtain symmetry between blocks and as shown by Taira and Colonius Taira_JCP_2007 (). The resulting matrix can be decomposed in the same way as performed by Perot Perot_JCP_1993 () using a th-order temporal approximation of the inverse of , i.e.

 A−1≈BN=N∑i=1Δti2i−1(M−1L)i−1M−1, (15)

where and . The approximation of the inverse has truncation error. Next, we employ the block-LU decomposition and arrive with three steps to solve the problem, i.e.

 Aq∗ =rn+bc1, (16) (GTBNGGTBNETnEnBNGEnBNETn)(ϕn+1~fn+1/2) =(GTEn)q∗−(−bc2uns), (17) qn+1=q∗ −BN(Gϕn+1+ETn~fn+1/2), (18)

where equation (16) is the so-called prediction step ( is the intermediate velocity flux), equation (17) is the pressure-Poisson step, and equation (18) is the projection step. The matrix in equation (16) and the block matrix on the left-hand-side of equation (17) are positive definite; hence the conjugate gradient method or the Cholesky factorization can be used to solve the linear systems in an efficient manner. Note that a modified finite volume scheme near the boundary may make the Laplacian non-symmetric. One may symmetrize by a similarity transform.

## 4 Numerical treatment of rigid body motion in fluid

### 4.1 Newton’s equations of motion

When the rigid body dynamics are not known a priori, it must be solved using Newton’s equations of motion. These equations in dimensionless form are repeated for convenience

 dusdt =1ρVs∮S¯¯¯¯¯¯¯τ⋅^ndS+(1−1ρ)g^eg, (19) d(¯¯¯¯Is⋅ωs)dt (20)

The Newton’s equations of motion in the IB framework pose a difficulty, because there is a fluid inside the solid body. One may consider the solid body as a region with fictitious boundary , which encompasses the IB forcing and the fluid within. Stresses over that given surface can be related to the flow field inside the volume and the volume forcing by

 ∮S+¯¯¯¯¯¯¯τ⋅^ndS =−∫V+~FdV+ddt∫V+udV, (21) ∮S+r×(¯¯¯¯¯¯¯τ⋅^n)dS =−∫V+~L×~FdV+ddt∫V+r×udV, (22)

where is the IB volume force density (applied on Eulerian grid), and the integration is over the fluid volume encompassed by surface . We note that the force density is non-zero only in a thin region close to the surface of the body. We assume that the spatial support of the discrete delta function is compact and consider the spatially limiting case as and , where and is the surface and volume of actual solid body, respectively. Following the suggestion by Kempe and Fröhlich kempe2012improved (), we evaluate the linear and angular acceleration of the fluid inside the particle as separate terms. A simpler approach suggested by Uhlman Uhlmann_JCP_2005 () is analyzed in section 7. Inserting equations (2122) in (1920), we arrive with Newton’s equations of motion in the IB framework

 ρVsdusdt =−∮SFdS+ddt∫VudV+Vs(ρ−1)g^eg, (23) ρd(¯¯¯¯Is⋅ωs)dt =−∮S~L×FdS+ddt∫Vr×udV, (24)

where we have substituted the volume force on the Eulerian grid with surface force on the Lagrangian grid by integrating over the delta function. Note that the first term on the right-hand side of equation (23) and (24) is integral over the surface of the solid body, which is later discretized using Lagrangian points, whereas the second term is volume integral over the fluid volume, which is later discretized using Eulerian mesh. We believe this to be the best representation of the equations in our method, since the unknown forcing is defined on the Lagrangian points, while unknown flow field is defined on the Eulerian mesh, as described is section 3.

While we only explain the detailed structure of matrices in two dimensions, the extension to three dimensions is straightforward. The moment of inertia tensor in a two-dimensional setting simplifies to a scalar, time independent constant . We introduce the solid body velocity variables , (translation in and directions, respectively) and (angular velocity) of center of mass for the rigid body. We can then discretize the Newton’s equations of motion (2324) as

 ρVsun+1s−unsΔt= −∑i∈L~fxi+dQn+1/2x+Vs(ρ−1)gx, (25) ρVsvn+1s−vnsΔt= −∑i∈L~fyi+dQn+1/2y+Vs(ρ−1)gy, (26) ρIsωn+1s−ωnsΔt= −∑i∈L(~Lnxi~fyi−~Lnyi~fxi)+dQn+1/2ω, (27)

where and are the relative coordinates of the body surface point measured from position of the center of mass , and and are the IB forcing components in and directions on the surface point , where . Here, is number of Lagrangian surface points. After the solution for the body velocity is found, the coordinates of Lagrangian points are updated as

 Ln+1xi= Lnxi+Δt(un+1s−ωn+1s~Lnyi), (28) Ln+1yi= Lnyi+Δt(vn+1s+ωn+1s~Lnxi), (29)

for all values from to . Note that we have omitted the time level index for the forcing, since it can have different indices depending on the particular integration scheme for coupling the equations of motion with the fluid equations (see next section). The terms and are the gravitational acceleration in and directions, respectively. We have used a first-order linear approximation of time derivative for volume integral of fluid inside particle. This derivative is denoted by , , for integrals , and , where and are the flow velocity components in and directions, respectively. The time level of the derivative is ; and the derivative depends on two previous flow fields and . The coefficients in front of integral value at each time level for first and higher-order approximations can be derived as explained by Fornberg fornberg1988generation (). For example, the time derivative approximation can be expanded as

 dQn+1/2x(un,un−1)=Qxun−Qxun−1Δt, (30)

where is the second-order midpoint rule for velocity component. In some fluid cells, the solid particle fills only a fraction of the cell. Therefore the quadrature in two dimensions takes the form

 Qxun=∑ijunijαijΔxΔy, (31)

where is the solid volume fraction in the fluid cell centered around grid point for fluid velocity component . Apart from solid volume fraction , equation (31) is the standard second-order midpoint rule ferziger2002computational (). As employed by Kempe and Fröhlich kempe2012improved (), the solid volume fraction can be found by the Heaviside step function and signed distance function between a grid point and the surface of the body as

 αij=∑4n=1−ζnH(−ζn)∑4n=1|ζn|, (32)

where summation is performed over all corners of the fluid cell . The term is assumed to be positive outside the body, and negative inside the body. We note that the Heaviside step function is used only as indicator function and is not discretized, contrary to delta function of IB formulation, which is discretized and smoothed.

In order to write the Newton’s equations of motion in a more compact form, we introduce a solid body velocity vector,

 (33)

a diagonal matrix of inertia222This is similar to generalized inertia matrix in the work of Wang and Eldredge wang2015strongly ().,

 IB=ρΔt⎛⎜⎝Vs000Vs000Is⎞⎟⎠, (34)

an inner-fluid integral vector,

 dQn+1/2B=(dQn+1/2xdQn+1/2ydQn+1/2ω)T, (35)

and a gravitational acceleration vector,

 GB=(gxgy0)T. (36)

We use the following arrangement of the IB forcing values,

 ~f=(~fx1…~fxnb−1~fxnb~fy1…~fynb−1~fynb)T. (37)

The force summation matrix333This is equivalent to transpose of rigid-body distribution operator  wang2015strongly (). at time level is defined as

 (38)

Using the above matrices and vectors, we can express the discrete Newton’s equations (2527) in the following algebraic form

 IB(un+1B−unB)=Nn+1/2B~f+dQn+1/2B+GB. (39)

Below we describe the explicit and implicit approaches for solving this equation in the framework of the IB projection method.

### 4.2 Explicit coupling between fluid and rigid-body dynamics

We start with describing forward Euler’s and Heun’s (predictor-corrector) methods for coupling the fluid-structure interaction, in order to compare with fully implicit scheme described in section 4.3. We denote the fluid solver as a function , which takes the solid body velocity vector and the coordinate vector as function arguments at time level , and returns the IB forcing . Note that we are reintroducing a specific time level for the IB forcing to illustrate the time level used in the coupling. Similarly, we denote the discrete Newton’s equations of motion (39) as a function , and coordinate update (2829) procedure as . Then, the Euler forward integration in time for the coupled system can be schematically written as

 ~fn+1/2 =NS(unB,Ln) (40) un+1B =ND(~fn), (41) Ln+1 =NC(unB), (42)

where is an interpolated IB forcing at time level . The outlined time stepping follows closely the definition of the forward Euler’s method, i.e. on the right hand side we only have parameters at time level .

Next, we describe Heun’s method (a predictor-corrector method). For the fluid-structure interaction problem, the first step of the Heun’s method is the forward Euler’s method discussed before

 ~f∗ =NS(unB,Ln), (43) u∗B =ND(~fn∗), (44) L∗ =NC(unB), (45)

where and are predictions for the solid body velocity and the coordinates respectively, and is the IB forcing at time level using the predicted value for interpolation. Next, the correction step is carried out for the Navier-Stokes equations and the results are used to advance the solid body dynamics

 ~fn+1/2 =NS(unB+u∗B2,Ln+L∗2), (46) un+1B =ND(~fn), (47) Ln+1 =NC(unB), (48)

where the IB forcing at time level is interpolated using and . Note that other combinations of time coupling are possible. For example, one could use the intermediate result from Navier-Stokes equations in Newton’s equations of motion or vice versa. We adopt the above approach, which resembles the Heun’s method as close as possible.

### 4.3 Implicit coupling between fluid and rigid body dynamics

We formulate implicit coupling (SC) as

 ~fn+1/2 =NS(un+1B,Ln) (49) un+1B =ND(~fn+1/2), (50) Ln+1 =NC(un+1B), (51)

where we observe that the output from solver (49) is used in (50) and vice versa. Therefore both the Navier-Stokes equations and the Newton’s equations of motion have to be solved simultaneously. Note that our SC scheme includes only implicit coupling for solid body velocity . There is no straightforward way to include body coordinates implicitly, because the interpolation and spreading operators and are depending on the coordinates of Lagrangian points , which would make the overall system non-linear. Consequently the overall temporal accuracy of the presented method is first order due to time-lagged interpolation and spreading operations444The dominant first-order contribution is .. In addition, since the system of equations now is solved at the same time, the force vector does not necessarily have to have a time level – it can be viewed as Lagrange multiplier. Nevertheless, we choose to keep the time level notation to be consistent with explicit methods.

The coupled Navier-Stokes and Newton’s equations can be written in the following algebraic form

 ⎛⎜ ⎜ ⎜ ⎜⎝A0GETn0IB0NnBGT000En(NnB)T00⎞⎟ ⎟ ⎟ ⎟⎠⎛⎜ ⎜ ⎜ ⎜ ⎜⎝qn+1un+1Bϕn+1/2~fn+1/2⎞⎟ ⎟ ⎟ ⎟ ⎟⎠=⎛⎜ ⎜ ⎜ ⎜⎝rnrnB0Δun+1B⎞⎟ ⎟ ⎟ ⎟⎠+⎛⎜ ⎜ ⎜⎝bc10−bc20⎞⎟ ⎟ ⎟⎠. (52)

Here the Newton’s equations of motion (39) are the second block-row and contains all known terms of the equations of rigid-body motion, i.e. the velocity at previous time step, the external forces at time level , and the derivative of velocity field inside the body at time level . The coupling between rigid-body dynamics and fluid dynamics is ensured using IB forcing both as volume forcing in Navier-Stokes equations (14) and as surface forcing for total force integral in Newton’s equations of motion (39). Furthermore, the coupling appears through the prescribed velocity at the solid body boundary (13), in which the velocity value at each boundary point is constructed using the unknown solid body velocity

 Enqn+1=−(NnB)Tun+1B+Δun+1B, (53)

where we have introduced as a prescribed difference between flow velocity at the boundary and the velocity of the body. For the no slip condition , whereas non-zero values correspond to a prescribed slip or penetration velocity, which results in some force exerted from the particle onto the fluid in order to match the prescribed slip or penetration velocity exactly. If the velocity difference is time dependent, the time level is selected to match the time level of solid body velocity.

Next, we can decompose the system (52) using a block-LU decomposition,

 (~A~Q~QT0)=(~A0~QT−~QT~BN~Q)(I~BN~Q0I), (54)

where we have defined

 ~A=(A00IB)    ~Q=(GETn0NnB). (55)

The matrix is symmetric and positive definite, since both and are symmetric and positive definite. The approximate inverse of order is

 ~BN=(BN00I−1B), (56)

where is the th-order approximation of from equation (15).

To summarize the method, let us list the three steps for IB projection method with implicit solver for rigid body dynamics;

1. [(I)]

2. prediction step

 (A00IB)(q∗u∗B)=(rn+bc1rnB), (57)
3. modified pressure Poisson solver

 (GTBNGGTBNETnEnBNGEnBNETn+(NnB)TI−1BNnB)(ϕn+1/2~fn+1/2) = (GTq∗+bc2Enq∗+(NnB)Tu∗B−Δun+1B), (58)
4. projection step (to enforce incompressibility and rigid body dynamics)

 (qn+1un+1B)=(q∗u∗B)−⎛⎝BNGϕn+1/2+BNETn~fn+1/2I−1BNnB~fn+1/2⎞⎠. (59)

It is noteworthy that the interaction between the solid body and the fluid is computed by modifying the Poisson matrix (58) using block matrix and by modifying the right-hand side using the predicted solid body velocity . Therefore, the added computational cost to the original method by Taira and Colonius Taira_JCP_2007 () is minimal. The size of modified pressure Poisson algebraic system is identical, while the algebraic system of prediction step is complemented only by or rows and columns for a single body in a two or three-dimensional setting, respectively (corresponds to a number of degrees of freedom for the rigid body).

It is straightforward to modify the method to include multiple bodies. In the case of bodies, Lagrangian points for all considered bodies must be assembled, and the corresponding velocity array, force array and diagonal matrix of inertia would be extended as

 ^uB=⎛⎜ ⎜ ⎜ ⎜ ⎜⎝u1Bu2B⋮umB⎞⎟ ⎟ ⎟ ⎟ ⎟⎠,    ^~f=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝~f1~f2⋮~fm⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,    ^IB=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝I1B0⋯00I2B⋯0⋮⋮⋱⋮00⋯ImB⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (60)

and the interpolation and force summation operators would be extended as

 ^ET=(ET1ET2⋯ETm),    ^NB=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝N1B0⋯00N2B⋯0⋮⋮⋱⋮00⋯NmB⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠, (61)

while keeping the structure of block-LU decomposition exactly the same. When describing dense particle suspensions one would need to complement the present method with appropriate collision model, such as the one employed by Kempe and Fröhlich kempe2012improved ().

## 5 Validation

In this section, all validation cases are computed using the implicit coupling, as described in the previous section. First, we carry out temporal and spatial convergence tests for a simple FSI problem. We then validate the present method on a freely falling/rising cylinder, and on migration of a neutrally buoyant cylinder in a shear flow.

### 5.1 Convergence

The convergence properties of the fluid solver and the IB method herein are reported by Perot Perot_JCP_1993 () and Taira and Colonius Taira_JCP_2007 (), respectively. We focus on the convergence when the fluid solver is coupled to rigid body dynamics. We select a simple two-dimensional problem – a circular cylinder with diameter falling under the influence of gravity in a fluid with kinematic viscosity (see Fig. 2a). The non-dimensional density of the cylinder is , the Reynolds number is , where is the terminal velocity, and the Gallileo number is .

To investigate the temporal convergence, we place the circular cylinder in the center of a square box with a width and a height of . We use a uniform mesh in the whole domain. For reference, we use a simulation with a moderate mesh spacing () and a small time step (). We carry out the simulation until , during which time the cylinder reaches the velocity . We perform a set of simulations for a range of time steps on the same mesh. We define an infinity norm as the maximum of error , i.e. . The infinity norm for the temporal error is shown in Fig. 2b, where it is observed that our method has the same convergence rate as the original method for practical time steps, i.e. third-order in time (in current test ), when a third-order approximation of the inverse of the Laplacian is used. In theory, the method is first-order in time due to the time-lagged interpolation and spreading operations. The first-order error starts to appear for smaller time steps (in current test , see Fig. 2b).

We found that the second-order approximation leads to unstable simulations for time steps larger than , which is similar to findings of Taira and Colonius Taira_JCP_2007 (); they showed that second-order approximation lacks positive-definiteness for higher time steps. We observed that the order of time derivative of integrals and plays a role in obtained convergence rate only for smaller time steps.

To investigate the spatial convergence, we reduce the box size to . For reference, we use a simulation with a very fine mesh () and a very small time step (). We perform the reference simulation until , during which the cylinder reaches the velocity . We use a range of mesh spacings , while keeping the time step constant. In order to compute the error , we interpolate the reference solution on the coarser grid using third-order interpolation in space. We define a L2-norm as , where is the number of points in the error field. As shown in Fig. 2c, the convergence rate is around in the L2-norm and around in the infinity norm, as observed by Taira and Colonius Taira_JCP_2007 () in their one-dimensional test case.

### 5.2 Freely falling and rising cylinder

There are a number of studies available for freely falling and rising bodies. For a review, we refer the reader to the work by Ern et al. ern2012wake (). For physically relevant range of Reynolds numbers Re (with density ratio close to unity), the falling cylinder problem has been investigated numerically by Namkoong et al. Namkoong_JFM_2008 (). They used an implicit coupling approach within a finite element method and adaptive body-fitted mesh with refined resolution in the cylinder wake.

For validation purposes, a density ratio and Reynolds number was selected. We carried out a simulation in a domain of size using a uniform mesh () with linearly expanding mesh at the boundaries, and a CFL condition . The boundary condition at the exterior of the simulation domain is no-slip. The vertical velocity of the falling cylinder is compared to the results from Namkoong_JFM_2008 () in Fig. 3a. The agreement is satisfactory, despite the relative simplicity of the current simulation method. The difference in the transient regime, where the wake instability develops, can be explained by the difference of rates at which numerical error accumulates and breaks the symmetry of cylinder wake. The transverse velocity of the falling cylinder is shown in Fig. 3b, and the trajectory of the falling cylinder is shown in Fig. 3c.

We present the vorticity field for freely falling and rising cylinders in Fig. 4 ( to with step ). We observe that the flow field of falling (a and b) cylinder is symmetric to that around a rising (c and d) cylinder. This is expected, since the Reynolds number for both configurations is the same and dimensionless densities ( and ) are very close to each other. Initially (), we see a symmetric vortex pair forming behind an accelerating cylinder. After the transient, flow perturbations have grown and the unstable symmetric solution has transitioned to the stable periodic vortex shedding state as seen in later times (). For the freely falling cylinder, we compare the Strouhal number – where is frequency of lift force oscillations or vortex shedding –, drag coefficient value – where is the drag force – and amplitude of lift coefficient oscillations with values reported by Namkoong et al. Namkoong_JFM_2008 () in Tab. 1, where we again observe satisfactory agreement. Although the flow fields are very similar, the freely rising cylinder has slightly higher Strouhal number, caused by the difference of the cylinder densities (Tab. 1); same observation was reported in Namkoong_JFM_2008 ().

### 5.3 Neutrally buoyant cylinder in shear flow

In order to further validate the present numerical method, we consider the problem of neutrally buoyant circular cylinder of diameter migrating in a Couette flow (see Fig. 5a). The channel height is with a uniform mesh in the direction. The channel width is with coordinate out of which the uniform mesh is used for (). The uniform mesh spacing is , and CFL number is set to . The upper and lower walls move in the -direction with velocities and , respectively, which gives a shear rate . The Reynolds number is , same as in the work of Feng et al. feng1994direct (), which we use as a reference. They used a body-fitted mesh with a finite element-solver. Their configuration is also used by Feng and Michaelides feng2004immersed (), Niu et al. niu2006momentum (), and Bhalla et al. Bhalla2013446 () for testing the coupling between fluid and rigid body dynamics.

We position the circular cylinder in the lower part of the channel () and release it with zero velocity and zero rotation. As observed in previous works feng1994direct (); feng2004immersed (); niu2006momentum (); Bhalla2013446 (), the cylinder rotates and migrates from the release position to the center of the channel (see Fig. 5b). The small oscillations in the trajectory (Fig. 5b) are due to the fact that Lagrangian force in IB methods is dependent on the position of Lagrangian points relative to the fluid grid. Breugem Breugem_JCP_2012 () has noted a similar behavior observed in IB simulations and referred to it as “grid locking”. Although we use the smooth 3-cell discrete delta function, the vertical migration is very slow compared to the rotation rate; and the error in vertical position associated to the “grid locking” can be seen. Consequently these oscillations illustrate the amplitude of error we have in the simulation.

We observe that rotation of the cylinder rapidly reaches a constant value, which is of the shear rate (see Fig. 6c), as reported in prior studies feng1994direct (); feng2004immersed (); niu2006momentum (); Bhalla2013446 (). To illustrate the flow field around the cylinder, we show equally spaced contour-lines of horizontal velocity from to (Fig. 6a,b). In case of a pure shear flow, the picture would consist of parallel lines only. When the cylinder is placed in the shear flow, initially – cylinder has not reached the rotation defined by the shear rate – there are significant distortions in the flow field (Fig. 6a). When the cylinder has reached the terminal rotation rate, the modifications of flow field are minor; the flow profile is compressed at the sides of the cylinder and flattened at the front and the back of the cylinder, which can be seen in Fig. 6b.

## 6 Numerical stability of coupling between rigid body dynamics and fluid dynamics

Let us examine the numerical stability of the current implicit method and compare it to an explicit coupling (see section 4.2) between a rigid body and the fluid. As a model problem, let us consider a body placed in a uniform free stream, which can move in the cross stream direction and freely rotate (2 degrees of freedom). We test the numerical method using circular cylinder with and without a splitter-plate of length clamped at the back of the cylinder. For the cylinder alone (see Fig. 7a), we expect vortex induced vibrations of the cylinder, a classical fluid-structure interaction problem investigated thoroughly in the literature sarpkaya2004critical (); williamson2004vortex (); Gabbai2005575 (). For the cylinder with the splitter plate (see Fig. 7b), we expect a drift caused by an inverted pendulum like (IPL) instability in addition to the VIV. As explained by Lācis et al. Lacis2014driftSymmetry (), the body orientation, when the splitter plate is parallel to the incoming free stream is always an equilibrium solution to the fluid-structure interaction problem. However, when the plate is sufficiently long, this solution becomes unstable in a manner similar to how an inverted pendulum becomes unstable under gravity. When this instability is triggered, the body turns until it reaches a new equilibrium turn angle, and it steadily drifts in the direction, in which the splitter plate has turned. We choose based on the free stream velocity , the cylinder diameter and the kinematic viscosity of fluid . There is no spring and no damper associated with the body.

We place the cylinder in the center of the computational domain of size . The mesh is uniform in the region for the flow around cylinder and for the flow around the cylinder with the splitter plate. The uniform grid spacing is and the CFL number is . The initial conditions shown in Fig. 8a and b have been obtained by performing simulations, in which we constrain all degrees of freedom for the body (fix the body) and let the flow evolve around it for time units.

In order to determine whether Euler’s forward method (see section 4.2) for coupling the rigid body and the fluid is stable for a given density ratio , we carry out the simulation from (stationary body) until the transient behavior has decayed and there are a number of periodic oscillations visible. Fig. 9a shows the vertical velocity for stable simulation of cylinder with .

We observe that after a short transient, the cylinder has reached a periodic transverse oscillatory state. Here, the so called VIV phenomenon sarpkaya2004critical (); williamson2004vortex (); Gabbai2005575 () has been reproduced – vortex formation at the lower side of the cylinder is shown in Fig. 10.

When the density ratio is reduced to , we observe that after some time (approximately time units) the scheme is unstable. The instability is illustrated in Fig. 9b by the oscillation energy

 Eosc=v2s+koscy2, (62)

where is an effective spring constant, and is the position of the body. The oscillation energy is based on kinetic energy of the cylinder and potential energy of a virtual spring; we neglect damping and forcing terms (for a complete VIV model see williamson2004vortex ()). We set the effective spring constant to value , such that the oscillation energy is the same, when the transverse velocity has maximal and zero values in the steady oscillation regime (). The critical density ratio is defined as the lowest value, for which we found the simulation to be stable, i.e. .

We carry out stability tests also for the remaining coupling algorithms, i.e. Heun’s method (section 4.2) and implicit method (section 4.3), for the cylinder with and without the splitter plate. All stability tests are performed over time units using an inverse Laplacian matrix with approximation order to save computational time. A summary of the obtained critical density ratios is listed in Tab. 2. It is observed that does not depend on the time step (in line with findings of Borazjani20087587 (), where they report no influence of time step on the numerical stability of vortex induced vibrations) and approximation order .

We can observe from Tab. 2 that adding a splitter plate behind a cylinder increases the value of critical density ratio