Shang-Gui Cai
email: shanggui.cai@utc.fr

Abdellatif Ouahsine

Julien Favier

Yannick Hoarau

Moving immersed boundary method

Keywords: Immersed boundary method, Projection method, Fractional step method, Implicit scheme, Lagrange multiplier, Incompressible viscous flow

come   L'archive ouverte pluridisciplinaire

INTRODUCTION

The immersed boundary method (IBM) has emerged in recent years as an alternative to traditional body-conforming mesh method for simulating fluid flows over complex and moving objects. Through adopting an appropriate boundary force in fluid equations for the presence of immersed solid boundaries, the simulations can be performed on a very simple Cartesian mesh. This significantly eases complicated mesh generations and eliminates moving boundary related issues, such as mesh distortions and mesh interpolation errors due to deforming-mesh and re-meshing.

Since first introduced by Peskin [START_REF] Peskin | Flow patterns around heart valves: A numerical method[END_REF] for modeling blood flow through a beating heart, IBM has been extended to various applications in scientific and engineering fields. In the original method, the immersed elastic membrane is represented by a series of massless Lagrangian markers where the boundary force is evaluated by using constitutive laws. Discretized delta functions are employed as kernel functions for the data exchange between the two independent meshes of fluid and solid. The immersed finite element method (IFEM) [START_REF] Wang | Extended immersed boundary method using FEM and RKPM[END_REF][START_REF] Zhang | Immersed finite element method[END_REF][START_REF] Liu | Immersed finite element method and its applications to biological systems[END_REF] was later developed in finite element formulations for general structures that occupy finite volumes within the fluid domain.

Previous methods are well suited for deformable solids owing to their physical basis, but the constitutive laws are generally not well posed when solids reach the rigid limit. Beyer and LeVeque [START_REF] Beyer | Analysis of a one-dimensional model for the immersed boundary method[END_REF] provided a solution by using a spring to attach the solids to an equilibrium location with a restoring force. Goldstein et al. [START_REF] Goldstein | Modeling a no-slip flow boundary with an external force field[END_REF] and Saiki and Biringen [START_REF] Saiki | Numerical simulation of a cylinder in uniform flow: Application of a virtual boundary method[END_REF] also proposed a feedback forcing strategy to control the velocity near the objects, which behaves as a system of springs and dampers. Nevertheless, artificial constants are introduced, which are ad hoc and should be chosen large enough in order to accurately impose the no-slip boundary condition. However large value makes the system very stiff and results in instabilities. The time step is severely limited, leading to a CFL number several magnitude smaller than the usual one [START_REF] Goldstein | Modeling a no-slip flow boundary with an external force field[END_REF][START_REF] Fadlun | Combined immersed-boundary finite-difference methods for threedimensional complex flow simulations[END_REF]. Mohd-Yosuf [START_REF] Mohd-Yosuf | Combined immersed Boundary/B-spline methods for simulation of flow in complex geometries[END_REF] and Fadlun et al. [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for threedimensional complex flow simulations[END_REF] proposed the direct forcing immersed boundary method to avoid the use of artificial constants via modifying the discrete momentum equation. No additional constraints are introduced to the time step. Instead of using the discrete delta function for velocity interpolation and force distribution, local velocity reconstruction approaches were employed to enforce the boundary condition. However Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] observed strong oscillations towards the boundary force. He attributed this problem to insufficient smoothing and re-used the discrete delta function in his direct forcing immersed boundary method. Although other strategies have also been proposed to enhance the local velocity reconstruction, special treatment should be taken for the phase change of cells near the moving boundaries [START_REF] Tseng | A ghost-cell immersed boundary method for flow in complex geometry[END_REF][START_REF] Wang | Algorithms for interface treatment and load computation in embedded boundary methods for fluid and fluid-structure interaction problems[END_REF][START_REF] Lakshminarayan | An embedded boundary framework for compressible turbulent flow and fluid-structure computations on structured and unstructured grids[END_REF].

In fact the boundary force is an unknown that is strongly coupled to the fluid velocity field. In the work of Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF], the boundary force is calculated explicitly by a tentative fluid velocity. The no-slip boundary condition can never be satisfied and large errors occur near the immersed boundaries. Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF] reduced the error by adding a forcing loop within a few iterations. Further improvement of the boundary condition imposition, however, requires numerous iterations for convergence as the multidirect forcing immersed boundary method [START_REF] Luo | Full-scale solutions to particle-laden flows: Multidirect Forcing and immersed boundary method[END_REF][START_REF] Breugem | A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows[END_REF]. Taira and Colonius [START_REF] Taira | The immersed boundary method: A projection approach[END_REF] proposed the implicit immersed boundary projection method (IBPM) by formulating the boundary force and the pressure into a modified Poisson equation and solving them simultaneously in an enlarged system with sophisticated solvers. Despite the mathematical completeness and rigour, IBPM may have convergence problems when an immersed boundary point is very close to a fluid grid point, as the singular property of the interpolation and distribution functions deteriorates significantly the condition number of the coefficient matrix of the original well-defined pressure Poisson equation (PPE) [START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF].

In this paper we propose the moving immersed boundary method (MIBM) to optimally maintain the accuracy of the implicit IBPM and the efficiency of the explicit direct forcing IBM. The projection method is served as the basic fluid solver where the proposed MIBM is integrated as a plug-in. Analogous to the role of the pressure in the projection method to satisfy the divergence-free condition, the boundary force is regarded as another Lagrange multiplier for the no-slip constraint in proposed MIBM. The global scheme follows the fractional step fashion and the fluid velocity, pressure and the boundary force are solved sequentially through the idea of operator splitting [START_REF] Cai | Improved implicit immersed boundary method via operator splitting[END_REF][START_REF] Cai | Implicit immersed boundary method for fluid-structure interaction[END_REF][START_REF] Cai | Computational fluid-structure interaction with the moving immersed boundary method[END_REF]. We follow the derivation of PPE in the projection method and derive an additional moving force equation for the boundary force. Therefore, the PPE is unchanged and immune from the convergence problem. Moreover, the force coefficient matrix is formulated to be symmetric and positive-definite so that generic linear system solvers can be applied directly.

The organization of this paper is as follows. First the fluid Navier-Stokes equations are discretized and a second order projection method is introduced as our fundamental fluid solver. In the following, the MIBM is presented in details and compared to other immersed boundary methods. In Section 5 a couple of numerical simulations are performed to validate the proposed MIBM. Finally, conclusions are drawn in Section 6.

FLUID EQUATIONS AND DISCRETIZATION

Consider the non-dimensionalized Navier-Stokes equations for an incompressible viscous fluid

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂u ∂t + ∇ • (u ⊗ u) = -∇p + 1 Re ∇ 2 u in Ω × [0, T ], ∇ • u = 0 in Ω × [0, T ],
u| Γ = w in [0, T ],

u| t=0 = u 0 in Ω, (1)
where u, p are the fluid velocity vector and the pressure. Re = UL/ν designates the Reynolds number, based on the reference velocity U , the reference length L and the kinematic viscosity ν. Directly solving above equations is very difficult. First the equations are non-linear due to the convective terms; Secondly there is no equation to compute the pressure directly; Moreover the pressure and the velocity are coupled through the continuity (incompressibility or divergence-free) condition, and the pressure is often regarded as a Lagrange multiplier to satisfy this constraint; Besides the solution of the pressure is not unique and is determined up to an additive constant. Numerical solutions will be discussed in this paper for overcoming these difficulties.

The above equations are discretized in space on a staggered mesh in order to prevent the socalled even/odd decoupling or checkerboard effect, as shown in Figure 1. The spatial derivatives are approximated by second-order central differences.

To discretize the equations in time, fully implicit scheme is superior to explicit one in terms of stability, which in turn requires non-linear iterations. This could be expensive in computation and the convergence is not always ensured. Non-linear iterations can be avoided by linearizing the convective terms, resulting in a non-symmetric coefficient matrix for the velocity. The matrix needs to be re-computed at each time step, which becomes very costly when the grid number increases. Fully explicit formulation seems to be very efficient as no iterations are needed. But the time step should be kept small enough to maintain stability. In two dimensions the constraints on the time step are the diffusive stability condition ∆t Re 2

1 ∆x 2 min + 1 ∆y 2 min -1 , (2)
and the convective stability condition of the usual CFL (Courant-Friedrichs-Lewy) type

∆t min ∆x min u max , ∆y min v max . (3)
It is easy to see that the diffusive constraint is more severe. Reducing the mesh size by half requires a four times smaller time step and it becomes more severe as the dimension increases. At low Reynolds number regime, the time constraint due to (2) dominates (3). It might be thought that for moderate to high Reynolds number flows, the diffusive stability condition is less restrictive. However in practice, the grid spacing is usually kept small under these circumstances for capturing small turbulence and the time step constraint of ( 2) is proportional to the square of the minimal mesh size.

In the present work, a semi-implicit time discretization scheme is employed, namely the convective terms are treated explicitly for avoiding the nonlinearity while the diffusive terms are treated implicitly for circumventing the severe diffusive time constraint. As a result, the entire system is linear and stable under the standard CFL condition. The velocity coefficient matrix remains symmetric and constant. To obtain a second order accurate system, we employ the second order Adams-Bashforth (AB2) scheme for the non-linear terms and the Crank-Nicolson (CN) scheme for the linear terms. The system now can be written as

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ u n+1 -u n ∆t + 3 2 N (u n ) - 1 2 N (u n-1 ) = -Gp n+1 + 1 2Re L(u n+1 + u n ), Du n+1 = 0, u n+1 | Γ = w n+1 , (4)
where L, N , G, D, are the discretized linear, non-linear, gradient, divergence operators, respectively. The superscript n + 1 and n represent the current time level and the past time level. The initial condition is hereafter omitted for convenience. The projection method, also refereed to fractional step method or time-splitting method, emerged in late 1960s as an effective tool to solve the pressure-velocity coupling problem, by splitting the system into a serial decoupled elliptic equations. The projection method is rooted in the Helmholtz-Hodge decomposition, which states that any smooth vector field v could be decomposed into the sum of a divergence-free part and a gradient of a potential field

v = v d + Gφ, (5)
where φ is often related to the pressure in the projection method. By taking the divergence of ( 5) and applying Dv d = 0, φ is the solution of the following Poisson equation

Lφ = Dv. ( 6
)
Once φ is calculated, the solenoidal velocity can be recovered by

v d = v -Gφ. (7)

Previous projection methods

The original projection method proposed by Chorin [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] and Témam [START_REF] Témam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (II)[END_REF] decouples the dynamic momentum equation from the kinematic incompressibility constraint by first estimating a tentative velocity û regardless of the pressure term, and then using the pressure to project the predicted velocity û into its solenoidal part u n+1 . The two sub-steps of prediction and projection are performed as

⎧ ⎪ ⎨ ⎪ ⎩ û -u n ∆t + 3 2 N (u n ) - 1 2 N (u n-1 ) = 1 2Re L(û + u n ), û| Γ = w n+1 , (8)
⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ u n+1 - û ∆t = -Gφ n+1 , Du n+1 = 0, u n+1 • n| Γ = w n+1 • n. (9)
By taking divergence of the first equation of ( 9) along with the incompressibility constraint, the actual realization of the projection step follows

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Lφ n+1 = 1 ∆t Dû, ∂φ n+1 ∂n | Γ = 0, u n+1 = û -∆tGφ n+1 , (10)
which is the same as ( 6) and [START_REF] Saiki | Numerical simulation of a cylinder in uniform flow: Application of a virtual boundary method[END_REF] when ∆t is absorbed to φ n+1 . In the projection method of Chorin [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] and Témam [START_REF] Témam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (II)[END_REF], the final pressure is set to p n+1 = φ n+1 . In spite of its efficiency, the original projection method suffers an irreducible splitting error of O(∆t), which deteriorates the original second order time discretization and prevents its extension to a higher-order method [START_REF] Perot | An Analysis of the Fractional Step Method[END_REF][START_REF] Liu | Projection method I: Convergence and numerical boundary layers[END_REF][START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. The error term can be found by adding the two sub-steps [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for threedimensional complex flow simulations[END_REF] and [START_REF] Mohd-Yosuf | Combined immersed Boundary/B-spline methods for simulation of flow in complex geometries[END_REF], and then comparing to (4)

1 2Re L(û -u n+1 ) = ∆t 2Re LGp n+1 , (11)
which is due to the time splitting scheme with the implicit treatment of the diffusive terms. Explicit treatment, however, would result in a severe limitation on the time step. It is rather natural to apply the physical boundary condition to the intermediate velocity û in the prediction step [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for threedimensional complex flow simulations[END_REF]. As a result, an artificial Neumann boundary condition ∂p n+1 /∂n| Γ = 0 is enforced on the pressure. This artificial homogeneous Neumann boundary condition introduces a numerical boundary layer to the solution, which prevents the method to be fully first-order [START_REF] Liu | Projection method I: Convergence and numerical boundary layers[END_REF][START_REF] Guermond | On stability and convergence of projection methods based on pressure Poisson equation[END_REF][START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. Improvements to the original projection method have been proposed in [START_REF] Goda | A multistep technique with implicit difference schemes for calculating two-or three-dimensional cavity flows[END_REF][START_REF] Braza | Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder[END_REF][START_REF] Van Kan | A second-order accurate pressure-correction scheme for viscous incompressible flow[END_REF] to achieve a higher order time accuracy by using an incremental scheme, which can be summarized as

⎧ ⎪ ⎨ ⎪ ⎩ û -u n ∆t + 3 2 N (u n ) - 1 2 N (u n-1 ) = 1 2Re L(û + u n ) -Gp n , û| Γ = w n+1 , (12)
⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ u n+1 - û ∆t = -Gφ n+1 , Du n+1 = 0, u n+1 • n| Γ = w n+1 • n, (13)
where an old value of pressure is retained in the prediction step, and φ n+1 here represents the pseudo pressure. The second sub-step is often performed as

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Lφ n+1 = 1 ∆t Dû, ∂φ n+1 ∂n | Γ = 0, u n+1 = û -∆tGφ n+1 , (14)
and the final pressure is updated by

p n+1 = p n + φ n+1 . ( 15
)
To study the spitting error, we sum up ( 12) and ( 13) and compare to [START_REF] Liu | Immersed finite element method and its applications to biological systems[END_REF]. Considering that the pseudo pressure is the approximation of φ n+1 = p n+1p n = ∆tp t , the splitting error is found to be of second order [START_REF] Perot | An Analysis of the Fractional Step Method[END_REF][START_REF] Armfield | An analysis and comparison of the time accuracy of fractional-step methods for the navierstokes equations on staggered grids[END_REF]

1 2Re L(û -u n+1 ) = ∆t 2Re LGφ n+1 = ∆t 2 2Re
LGp t .

Note that the physical boundary condition is still assigned to the intermediate velocity in the prediction step [START_REF] Wang | Algorithms for interface treatment and load computation in embedded boundary methods for fluid and fluid-structure interaction problems[END_REF]

∂n | Γ = ∂p n ∂n | Γ = • • • = ∂p 0 ∂n | Γ , (17)
is enforced on the final pressure. This pressure boundary condition is not physical, thus it introduces a numerical boundary layer and prevents the scheme to be fully second order [START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. This error is irreducible, hence using a higher order time stepping scheme will not improve the overall accuracy.

Rotational incremental pressure-correction projection method

To obtain a solution of second order accuracy with consistent boundary conditions, we propose to use the rotational incremental pressure-correction projection method of [START_REF] Timmermans | An approximate projection scheme for incompressible flow using spectral elements[END_REF][START_REF] Guermond | An overview of projection methods for incompressible flows[END_REF]. The essential idea of this method is to absorb the splitting error into the pressure so that the sum of the substeps is consistent with the original discretized momentum equation [START_REF] Liu | Immersed finite element method and its applications to biological systems[END_REF]. By considering the identity ∇ 2 u = ∇(∇ • u) -∇ × ∇ × u, the error term ( 16) can be rewritten as

1 2Re L(û -u n+1 ) = 1 2Re G(Dû), (18)
where ∇ × ∇ × û = ∇ × ∇ × u n+1 is used, which can be verified by the Helmholtz-Hodge decomposition. Now the error term in this form can be absorbed into the pressure

p n+1 = p n + φ n+1 - 1 2Re
Dû.

Most importantly, the pressure boundary condition is consistent with the original system. Therefore, no numerical boundary layer will be generated with this scheme. Higher than second order accuracy can be achieved if a higher-order time-stepping scheme is used. In the present work, the second order accuracy is found to be sufficient with the AB2 scheme and the CN scheme. The overall rotational incremental pressure-correction projection method can be summarized as follows

⎧ ⎪ ⎨ ⎪ ⎩ û -u n ∆t + 3 2 N (u n ) - 1 2 N (u n-1 ) = 1 2Re L(û + u n ) -Gp n , û| Γ = w n+1 , (20)
⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Lφ n+1 = 1 ∆t Dû, ∂φ n+1 ∂n | Γ = 0, u n+1 = û -∆tGφ n+1 , (21)
p n+1 = p n + φ n+1 - 1 2Re
Dû.

Pressure Poisson equation solver and Parallel computing

The aforementioned discretized equations lead to a set of linear systems to be solved. Among them the pressure Poisson equation is the most time-consuming part due to its high condition number, which is generally solved iteratively to save computational time and storage. (Bi-CGSTAB), generalized minimum residual (GMRES), etc., are very efficient for this problem.

In addition, more efficiency can be achieved if pre-conditioning is applied, such as the incomplete Cholesky (IC) factorization, the incomplete lower-upper (ILU) decomposition and the approximate inverse (AINV) [START_REF] Chow | Approximate inverse preconditioners via sparse-sparse iterations[END_REF]. The MG method is found to be more efficient when used as a pre-conditioner in conjunction with Krylov solvers instead of a pure solver.

To further improve this work, the code is extended to allow parallel computing. First we integrate our method into the PETSc library [START_REF] Balay | PETSc users manual[END_REF], which employs the MPI for communications between the CPU cores. In the second mode, we parallelize the code through using the CUDA CUSP library on GPU [START_REF] Dalton | Cusp: Generic parallel algorithms for sparse matrix and graph computations[END_REF]. In fact the CPU consists of a few cores optimized for sequential serial task, while the GPU may have massive smaller cores at the same price which is extremely efficient for handling multiple tasks simultaneously. Therefore, we send the parallelable and computationally intensive parts of the application to the GPU and run the remainders on the CPU. From the practical point of view, the second mode runs significantly fast. Table I illustrates the performances of the two parallelization modes. The test is performed by solving the PPE on a 400 × 400 grid with the Neumann boundary condition applied at all the boundaries. As a matter of fact, this system does not possess a unique solution (up to an additional constant). We pin a fixed value to one cell to remove the zero eigenvalue, as suggested in [START_REF] Taira | The immersed boundary method: A projection approach[END_REF][START_REF] Kassiotis | Nonlinear fluid-structure interaction problem. Part I: implicit partitioned algorithm, nonlinear stability proof and validation examples[END_REF]. The calculation is done on the platform PILCAM2 with the CPU Intel Xeon X7542 and the GPU Quadroplex 2200 S4. The process time decreases approximately by half when we double the cores of CPU from 1 to 16 in the test with CG solver. About 1.5-2.8 times' acceleration has been achieved when the MG method is applied as a preconditioner for the CPU parallelization. The parallelization of GPU greatly accelerates the calculation up to 40 times in the test with MG preconditioner. Different preconditioners like the AINV and the MG are compared in Table II

φ 4 (r) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1 8 3 -2|r| + 1 + 4|r| -4r 2 , |r| < 1, 1 8 5 -2|r| --7 + 12|r| -4r 2 , 1 |r| < 2, 0, otherwise, (30)
which is widely used in the literature. Roma et al. [START_REF] Roma | An adaptive version of the immersed boundary method[END_REF] also designed a 3-point-width function specially for the staggered mesh [START_REF] Beyer | Analysis of a one-dimensional model for the immersed boundary method[END_REF]; ----, the 3-point-width function of Roma et al. [START_REF] Roma | An adaptive version of the immersed boundary method[END_REF]; --, the 4-point-width function of Peskin [START_REF] Peskin | The immersed boundary method[END_REF].

φ 3 (r) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ 1 3 1 + -3r 2 + 1 , |r| < 0.5, 1 6 5 -3|r| --3(1 -|r|) 2 + 1 , 0.5 |r| < 1.5, 0, otherwise. ( 31
) -3 -2 -1 0 1 2 3 0 0.2 0.4 0.6 0.8 1 Figure 4. Comparison of the one-dimensional function φ(r). • • • •, the 2-point-width hat function
The functions are plotted and compared in Figure 4. The one-dimensional function of Roma et al. [START_REF] Roma | An adaptive version of the immersed boundary method[END_REF] has a relative smaller support than the four-point version of Peskin [START_REF] Peskin | The immersed boundary method[END_REF], providing a sharper interface and a better numerical efficiency while maintaining good smoothing properties. The discrete delta functions used in the present work have the following properties:

• δ h has a narrow support to reduce the computational cost and to obtain a better resolution of the immersed boundary.

• δ h is second order accurate for smooth fields.

• δ h satisfies certain moment conditions to meet the translation invariant interpolation rule, namely the total force and torque are equivalent between the Lagrangian and Eulerian locations.

x∈gh δ h (x -X)h 2 = 1 (zeroth moment condition), (32)
ũ = u * + ∆tf n+1 . (38)
(6) Implicit treatment of the viscous term

1 ∆t û - 1 2Re Lû = 1 ∆t ũ. (39)
(7) Project the fluid velocity into the divergence-free field and update the pressure

Lφ n+1 = 1 ∆t Dû, (40)
u n+1 = û -∆tGφ n+1 , (41)
p n+1 = p n + φ n+1 - 1 2Re
Dû.

Here u * , ũ, û, u n+1 represent the fluid velocity at each stage of the fractional step method, i.e., the prediction step of explicit terms, the immersed boundary forcing step, the viscous prediction step, and the projection step. U b (X l ) is the solid velocity of the lth element at the immersed boundary.

The method of Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] is favored in the literature as it is computational inexpensive due to its explicit treatment of the boundary force. However, numerical simulations have shown that it fails to impose the velocity boundary condition exactly on the immersed boundary [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF]. A forcing error is introduced which is irreducible and depends on the time step and Reynolds number Re. This error comes from the fact that the tentative fluid velocity u * is used for the boundary force evaluation. The ideal velocity should be the final fluid velocity u n+1 while it is unknown at the immersed boundary forcing step. Otherwise we need to iterate the whole system to achieve u n+1 implicitly, which could be too cumbersome to perform. But this implies one way of reducing the forcing error by choosing the closest value to the final velocity.

Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF] suggested to perform the viscous prediction step first and then use the intermediate velocity û to compute the boundary force. To further improve the accuracy, a forcing loop is added in the immersed boundary forcing step. This additional loop is performed within few iterations without convergence. The method of Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF] can be expressed as (1) Prediction of the explicit terms

u * = u n + ∆t - 3 2 N (u n ) - 1 2 N (u n-1 ) -Gp n + 1 2Re Lu n . (43)
(2) Viscous prediction step

1 ∆t û - 1 2Re Lû = 1 ∆t u * . (44)
(3) Immersed boundary forcing loop Loop for k = 1 to 3 with

û(0) = û Û(k) (X l ) = nx i=1 ny j=1 û(k-1) δ h (x i,j -X l )h 2 , (45)
(k) (X l ) = U n+1 b (X l ) -Û(k) (X l ) ∆t , (46)
f (k) (x i,j ) = nb l=1 F (k) (X l )δ h (x i,j -X l )∆V l , (47)
ũ(k) = û(k) + ∆tf (k) , (48)
û(k) = ũ(k) . (49)
End loop (4) Projection step and update of the final fields

Lφ n+1 = 1 ∆t Dũ, (50)
u n+1 = ũ -∆tGφ n+1 , (51)
p n+1 = p n + φ n+1 - 1 2Re
Dû.

If full convergence of the forcing loop is required more iterations are needed, such as the multidirect forcing scheme of [START_REF] Luo | Full-scale solutions to particle-laden flows: Multidirect Forcing and immersed boundary method[END_REF][START_REF] Breugem | A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows[END_REF]. However, the convergence rate of this iteration becomes very slow after several iterations. The computational cost increases hugely when more Lagrangian points are involved in the additional forcing loop. Therefore the number of iteration is usually kept low for the computational efficiency. Even though the error is reduced, the method of Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF] is still explicit. The exact no-slip boundary condition can never be satisfied.

To impose the no-slip boundary condition exactly, Taira and Colonius [START_REF] Taira | The immersed boundary method: A projection approach[END_REF] proposed the implicit immersed boundary projection method (IBPM) by combining the boundary force and the pressure into a modified Poisson equation and solving them simultaneously in one single projection step. However convergence problem may occur as one boundary point is very close to a fluid grid point [START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF], because the singular property of the interpolation and distribution functions undermines the coefficient matrix condition number of the PPE. Ji et al. [START_REF] Ji | A novel iterative direct-forcing immersed boundary method and its finite volume applications[END_REF] proposed to iterate each part to get rid of the convergence problem, which is inevitable computational expensive.

Novel implicit immersed boundary method

In this subsection we present a novel implicit but efficient IBM variant, termed as the moving immersed boundary method (MIBM) in this paper. The objective of MIBM to maintain the efficiency of the explicit direct forcing IBM but with an improved accuracy like the multidirect forcing IBM and the IBPM.

To this end, we first take the immersed boundary forcing part from the explicit IBM of Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF] for consideration, i.e., [START_REF] Tritton | Experiments on the flow past a circular cylinder at low Reynolds numbers[END_REF], [START_REF] Wang | An immersed boundary method based on discrete stream function formulation for two-and three-dimensional incompressible flows[END_REF], [START_REF] Lai | An immersed boundary method with formal second-order accuracy and reduced numerical viscosity[END_REF] and [START_REF] Williamson | Oblique and parallel modes of vortex shedding in the wake of a circular cylinder at low Reynolds numbers[END_REF]. By dropping the superscripts for convenience, the immersed boundary forcing part is written as

Û = T û, (53)
F = U b - Û ∆t , (54)
ũ = û + ∆tf . ( (55)
) 56
We require that the interpolated velocity satisfies the no-slip wall boundary condition on the immersed interface after the immersed boundary forcing, namely T ũ = U b , then

T (û + ∆tf ) = U b . (57)
Substituting ( 55) into [START_REF] Wang | Two dimensional mechanism for insect hovering[END_REF] gives

T (û + ∆tSF) = U b , (58)
which can be rearranged in order to separate the boundary force

(T S)F = U b -T û ∆t . ( 59
)
We donate M = T S the moving force coefficient matrix. M is a function of the solid position, which changes its value as the boundary moves. Thus the force is redistributed just like the boundary force moves. The moving force equation can be rewritten in a more concise form

MF = F e , (60)
where F e = (U b -T û)/∆t is exactly the explicit forcing value used in [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF].

Compared to the modified Poisson equation in the IBPM of [START_REF] Taira | The immersed boundary method: A projection approach[END_REF], the moving force equation ( 60) is much smaller in size and easier to work with. At each dimension (x or y), the size of the force coefficient matrix is n b × n b since T ∈ R nb×nxny and S ∈ R nxny×nb . Therefore, for moving boundaries, its update is computational less expensive than the modified Poisson equation.

Note that S = (∆V l /h 2 )T T if the same function is used for interpolation and spreading, where ∆V l /h 2 ≈ 1 is the volume ratio between the fluid and the solid cell. As a result, the moving force coefficient matrix M = (∆V l /h 2 )T T T is symmetric. It is also found that M is positive-definite irrespective of the time step and the approximation order as in the IBPM [START_REF] Taira | The immersed boundary method: A projection approach[END_REF]. Moreover, the moving force equation is well conditioned, which converges quickly by using the conjugate gradient method. Now we incorporate this moving force equation into the rotational incremental pressurecorrection projection method. For the sake of simplicity, we rewrite the governing equations [START_REF] Liu | Projection method I: Convergence and numerical boundary layers[END_REF] as

u n+1 -u n ∆t = H + P + F, (61)
Du n+1 = 0, (62)
T u n+1 = U n+1 b , (63)
where H, P and F are the operators defined as

H := - 3 2 N (u n ) - 1 2 N (u n-1 ) + 1 2Re L(u n+1 + u n ) -Gp n , (64)
P := -Gφ n+1 , (65)
F := SF n+1 . ( 66
)
To decouple the momentum equation (61) from the divergence free condition (62) and the no-slip wall condition on the interface (63), we perform the following operator splitting algorithm:

(1) Prediction step by ignoring the immersed objects

û -u n ∆t = H(û). (67)
(2) Immersed boundary forcing step for satisfying the no-slip wall condition on the interface

ũ - û ∆t = F, (68)
T ũ = U n+1 b . (69)
Applying ( 69) to (68) gives the moving force equation that we have defined previously

MF n+1 = U n+1 b -T û ∆t . ( 70
)
Once the boundary force is determined, we correct the fluid velocity with

ũ = û + ∆tSF n+1 . (71)
(3) Projection step for obtaining the divergence free velocity u n+1 and the final pressure p n+1

u n+1 - ũ ∆t = P, (72)
Du n+1 = 0. ( 73
)
Applying the divergence operator to (72) and using the divergence free condition (73) gives

Lφ n+1 = 1 ∆t Dũ, (74)
u n+1 = ũ -∆tGφ n+1 . ( 75
)
The final pressure is advanced by

p n+1 = p n + φ n+1 - 1 2Re
Dû.

(76)

Figure 6 shows the global structure of MIBM. The overall scheme follows the regular fractional step method so that the velocity, the pressure and the force are decoupled. Even though the interface velocity condition is enforced before the projection step, we have found that the velocity on the immersed boundary is essentially unchanged after the projection step. The same observation has also been made by Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF] and Fadlun et al. [START_REF] Fadlun | Combined immersed-boundary finite-difference methods for threedimensional complex flow simulations[END_REF]. It is worth noting that the present MIBM recovers to the explicit method of Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF] with one iteration in the forcing loop, if M is set to the identity matrix. However it is not the case, hence our method is implicit.

Comparison of performance

To demonstrate the accuracy and efficiency of present moving immersed boundary method, we perform the following test Given u 0 (x, y) = e x cos y -2, 0 x, y 1,

Find F such that u(x, y) = u 0 (x, y) + ∆tSF = U b on Γ s ,
where Γ s is described with a circle of a radius of 0.2 at (0.52, 0.54) and U b = 0. The domain is covered by 64 × 64 nodes with around 81 Lagrangian points on the circle surface. ∆t is set to 1.

In this test, the fluid equations are not solved and only the immersed boundary forcing part is considered. The initial field u 0 (x, y) can be seen as a predicted fluid velocity component in one direction. This test is to examine different forcing strategies for imposing the desired velocity U b at the interface Γ s via a boundary force F . To facilitate the accuracy study, we define the velocity error norms of L 2 and L ∞ as follows for i = 1, . . . , n x , j = 1, . . . , n y where u ref represents the reference value. It is worth noticing that the L 2 -norm is a good measure of the global error while the L ∞ -norm provides a good indicator for the local error. Figure 7a displays the result of the explicit direct forcing IBM of Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF], where u is far away from zero over the immersed boundary compared to Figure 7c. The accuracy is improved after 3 iterations with the method of Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF], as shown in Figure 7b. Figure 7d reveals that the results are nearly the same for present MIBM with the iterative multidirect forcing IBM of Luo et al. [START_REF] Luo | Full-scale solutions to particle-laden flows: Multidirect Forcing and immersed boundary method[END_REF] and Breugem [START_REF] Breugem | A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows[END_REF]. Table III compares the computational time and velocity error on the interface of these immersed boundary methods. The error is measured in L 2 -norm and the tolerance is 1 × 10 -15 . The method of Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] is the quickest due to its explicit nature, but it suffers a large error of 3.01 × 10 -1 on the immersed interface. The forcing loop of Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF] reduces the error by a factor of 4 with 3 iterations. However, the error of 7.41 × 10 -2 is still considered large.

||e u || 2 = 1 n x n y nx i=1 ny j=1 (u i,j -u ref i,j ) 2 1/2 , ( 77
)
||e u || ∞ = max|u i,j -u ref i,j |, (78)
The iterative multiforcing IBM of Luo et al. [START_REF] Luo | Full-scale solutions to particle-laden flows: Multidirect Forcing and immersed boundary method[END_REF] and Breugem [START_REF] Breugem | A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows[END_REF] is required to converge towards the machine precision, but it takes approximately 606 times more additional computational effort than the explicit method of Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF]. Actually, the convergence rate in the multiforcing IBM decreases dramatically after about 10 iterations, as shown in Figure 8. In order to reduce the error to 1 × 10 -6 around 1000 iterations are needed and 4443 iterations for the machine precision. The present MIBM converges to the same machine precision only with 60 iterations by using the conjugate gradient solver. The iteration can be further reduced if preconditioning is taken, but we find that the conjugate gradient solver is sufficient for fast convergence. The computation is not increased considerably compared to the explicit method of Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF], as we can see that the present method only takes twice the amount of computational time of the direct forcing IBM of Uhlmann [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF]. It also worth noticing that present MIBM is almost as efficient as the method of Kempe and Fröhlich [START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF].

Iteration Velocity Error

Figure 8. Comparison of convergence between present MIBM (--) and the multidirect forcing IBM of Luo et al. [START_REF] Luo | Full-scale solutions to particle-laden flows: Multidirect Forcing and immersed boundary method[END_REF] and Breugem [START_REF] Breugem | A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows[END_REF] (----).

RESULTS

Taylor-Green vortices

We first consider the two-dimensional unsteady case of an array of decaying vortices to assess the accuracy of the fluid solver. The analytical solution of the Taylor-Green vortices is given by u(x, y, t) = -cos(πx)sin(πy)e -2π 2 t/Re , v(x, y, t) = sin(πx)cos(πy)e -2π This simulation is performed on a square domain Ω = [-1.5, 1.5] × [-1.5, 1.5] and the Reynolds number Re is prescribed to 10. The initial and boundary conditions are provided by the exact solution. We advance the equations for 0 t 0.2.

To study the temporal accuracy, we compare the results at t = 0.2 to a reference solution obtained by a very fine time step ∆t = 1 × 10 -4 with the spatial resolution of ∆x = ∆y = 9.375 × 10 -3 . The errors on the velocity component u are computed by subtracting the reference solution from other numerical solutions (∆t ∈ [0.00125, 0.01]), to cancel out the error due to spatial discretization. The L 2 , L ∞ error norms are then displayed in Figure 9a on a log-log plot. A second order temporal accuracy is observed, which confirms previous error estimation analysis for the rotational incremental pressure-correction projection method.

We also expect a second order spatial accuracy since the second order central differencing scheme is used for all the derivatives in this case. We use a small time step ∆t = 1 × 10 -4 to ensure that the temporal discretization error is negligible compared to the spatial one, and then vary the computational grids (n x × n y = 20 × 20, 40 × 40, 80 × 80, and 160 × 160). The error is obtained by comparing the results to the analytical solution. Figure 9b shows the spatial discretization error, indicating a second order spatial accuracy. It is well known that the discrete delta function undermines the space accuracy of the original fluid solver. Now we embed a circular cylinder of a unit radius in the center of the computational domain to study the accuracy of our MIBM. The time dependant no-slip boundary condition at the immersed cylinder surfaces is enforced by current MIBM. Figure 9b shows the variation of the velocity error as a function of the mesh size. It is evident that current MIBM introduces errors the original fluid solver but it still retains the second order accuracy, which corresponds to the interpolation properties of the discrete delta function for smooth fields.

Lid-driven cavity flow with an embedded cylinder

In this test, we compare current immersed boundary method with the traditional body-conforming mesh method. The domain configuration and the boundary conditions are taken the same as in the classical lid-driven cavity flow case, namely the top wall is moving with a constant velocity u ∞ = 1 while the others are stationary walls, except that we place a cylinder in the domain center. In order to compare with Vanella and Balaras [START_REF] Vanella | A moving-least-squares reconstruction for embedded-boundary formlations[END_REF], the diameter of the cylinder is set to D = 0.4L with L being the cavity length. The Reynolds number is 1000 in this study based on the cavity length. A uniform mesh of 200 × 200 is employed in the immersed boundary method, and the same mesh size is used for the body-conforming mesh method for comparison. Body-conforming mesh method (0.6906, 0.6872) (0.0791, 0.0721) (0.8849, 0.1063)

Table IV. Comparison of vortices center positions for the proposed immersed boundary method and the body-conforming mesh method, where (x 1 , y 1 ), (x 2 , y 2 ), (x 3 , y 3 ) are the vortices centers at the upper right to the cylinder, at the lower left corner and at the lower right corner respectively.

The flow reaches a final steady state as the time advances. Figure 10 shows the vorticity contours and streamlines for the flow at Re = 1000, which are similar to the results of [START_REF] Vanella | A moving-least-squares reconstruction for embedded-boundary formlations[END_REF]. As we can see, three vortices emerge in the flow. One at the upper right position of the cylinder and two near the bottom at each corners. It is noteworthy that the upper vortex is generated by the presence of the cylinder. The flow fields outside the cylinder are essentially the same for current MIBM and the body-conforming mesh method. The only difference is that there is a flow inside the cylinder in the immersed boundary method, which however is the key idea of the immersed boundary method to replace the solid domain with fluid. The velocity component u at the vertical midline x = 0.5 and the velocity component v at the horizontal midline y = 0.5 are plotted in Figure 11. The velocity profiles of both methods match pretty well. The location of the three vortices centers are also listed in Table IV. Very close results have been obtained. Next we study the grid convergence for assessing the accuracy of present method for nonsmoothed field. A series of computations are performed on a hierarchy of grids (70 × 70, 90 × 90, 126 × 126, 210 × 210 and 630 × 630). The variation of error of the velocity component u along with the grid spacing is displayed in Figure 12, showing a convergence rate of about 1.13. This is because the flow becomes not smooth near the immersed surface in this case, and the discrete [START_REF] Beyer | Analysis of a one-dimensional model for the immersed boundary method[END_REF] analysed various discrete delta functions and pointed out that the second order accuracy can be recovered through using different functions for interpolation and spreading. This results in non-symmetric coefficient matrix of the boundary force in MIBM, which can be solved with the GMRES or Bi-CGSTAB methods.

Flow over a stationary circular cylinder

The flow past a stationary circular cylinder is considered as a canonical test case to validate current method, since a great amount of experimental and numerical studies at different Reynolds numbers are available for comparison. The flow characteristics depend on the Reynolds number Re = u ∞ D/ν, based on the inflow velocity u ∞ , the cylinder diameter D = 1 and the fluid kinematic viscosity ν. The simulation is performed in a rectangular domain, where the fluid flows from the left to the right (see Figure 13). At left boundary, a uniform velocity of u ∞ = 1 is imposed; The free slip boundary conditions are applied at lateral boundaries; At outlet, the convective boundary condition ∂u/∂t + u ∞ ∂u/∂x = 0 is employed for reducing the reflection effects because of the finite artificially truncated domain. The cylinder is placed at the center of the computational domain. The fluid domain is covered with a uniform mesh, and the cylinder surface is represented by a set of uniformly distributed Lagrangian points with δs ≈ h.

For comparison the drag and lift coefficients are defined as

C D = F D 1 2 ρu 2 ∞ D , C L = F L 1 2 ρu 2 ∞ D , (80)
where F D , F L are the drag and lift forces on the cylinder exerted by the fluid, respectively. The fluid density ρ is set to 1 here. As a matter of fact, the spreading and interpolation operators constructed from the regularized delta function conserve the total force, hence F D and F L can be computed directly by summing up the forces over all the Lagrangian points

F D F L = - nb l=1 F(X l )∆V l . (81) D Inflow u ∞ Convective outlet
Free-slip Free-slip The time-averaged values of the wall vorticity ω z and the wall pressure coefficient C P are shown in Figure 21 for Re = 100. Good agreements have been found compared to the results of Braza et al. [START_REF] Braza | Numerical study and physical analysis of the pressure and velocity fields in the near wake of a circular cylinder[END_REF]. The effects of different discrete delta functions on the results are also tested in for Re = 100, 200, where a domain of Ω = 30D × 30D is used and the mesh resolution is set to h = 0.029D.

(Ω = 30D × 30D, h = 0.04D) 1.355 ±0.042 ±0.677 0.200 Present (Ω = 30D × 30D, h = 0.029D) 1.365 ±0.044 ±0.696 0.200 Present (Ω = 30D × 30D, h = 0.
A careful grid convergence study is also performed to examine the order of accuracy in this case. Since the exact solution does not exist, we use the solution calculated on a highly resolved grid of 630 × 630 as our reference for computing the error. The computation domain is taken as [-2D, 2D] × [-2D, 2D] with the Reynolds number Re = 100. The equations are advanced until 0.2 and a relative small time step of 5 × 10 -4 is chosen such that the time discretization error will not influence the results. Same computations but on different grids are performed and compared the reference solution, namely 45 × 45, 70 × 70, 90 × 90, 126 × 126 and 210 × 210. The distribution of velocity error in the x-direction for the 90 × 90 grid is shown in Figure 22. Large magnitudes of error in velocity are located near the cylinder. Figure 23 displays the L 2 norm of this error on a log-log plot. A convergence rate of around 1.21 is observed.

Re = 1000

We further extend our method to a higher Reynolds number flow Re = 1000. At this regime, the convection effects become predominant and the boundary layer thickness decreases, which can be estimated by δ ≈ D/ √ Re = 0.032. To capture the thin boundary layer, a fine grid resolution of h = 0.01D is taken, as recommended in [START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF][START_REF] Apte | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF]. Note that the grid resolution is only marginal for resolving the boundary layer at this Reynolds number. Nevertheless, the results are satisfactory and the essential features of the flow are well captured. The computational domain is chosen to be [-20D, 20D] × [-20D, 20D]. The two-point-width hat function φ 2 is employed in this case as it provides a sharp interface. Figure 25 shows the instantaneous vorticity field. The coefficients of drag and lift are plotted in Figure 24. Note that the flow is inherently three-dimensional at this Reynolds number. We compare our simulations with other two-dimensional results available in the literature. The properties of the drag and lift coefficients are summarized in Table VIII. Good agreements have been found. The computational domain is chosen to be 14D × 14D, as shown in Figure 26. The cylinder is initially located at the center of the computational domain. The outflow boundary condition ∂u/∂n = 0 is applied at the domain contours. A uniform mesh of 560 × 560 is adopted for the fluid domain and the cylinder is represented by 126 points due to δ s ≈ h. The transient no-slip velocity boundary condition at the cylinder surface is enforced by present MIBM at each time level

u(t) = -2πfA cos(2πft). (83)
The pressure and vorticity contours at four different phases (φ = 2πft = 0 • , 96 • , 192 • , 288 • ) are shown in Figure 27, where two counter-rotating vortices are formulated during the oscillation. The vortices contours are drawn from -3 to 3 with an increment of 0.4, which display the same structure as in [START_REF] Dütsch | Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan-Carpenter numbers[END_REF].

Figure 28 shows the profiles of the velocity components u and v at four different streamwise locations (x = -0.6D, 0D, 0.6D, 1.2D) for three phase (φ = 2πft = 180 • , 210 • , 330 • ). The experimental results of [START_REF] Dütsch | Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan-Carpenter numbers[END_REF] by LDA measurements are also plotted for comparison. The velocity profiles outside the cylinder agree well those of [START_REF] Dütsch | Low-Reynolds-number flow around an oscillating circular cylinder at low Keulegan-Carpenter numbers[END_REF]. The only discrepancy is the velocity inside the cylinder. Since the present IBM treats the solid domain as fluid, the velocity is non-zero inside the cylinder. From Figure 28 we can see that this treatment, however, does not influence the flow field outside the solid. Various internal treatments of the body have been discussed in the work of [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF], such as applying the force inside the body and thus changing the velocity distribution. Iaccarino and Verzicco [START_REF] Iaccarino | Immersed boundary technique for turbulent flow simulations[END_REF] also concluded that for direct forcing IBM, there is essentially no difference. Therefore, for simple implementation we just leave the interior of the solid free to develop a flow without imposing anything.

Flow around a flapping wing

In this example, we investigate the flow induced by a flapping wing, in order to demonstrate the ability of current method for handling non-circular object in both translational and rotational motions. The configuration of this problem is shown in Figure 30. The hovering wing is a geometrical 2D ellipse with major axis c (chord length) and minor axis b. The aspect ratio is defined as e = c/b. The wing is initially located at the origin with an angle of attack of θ 0 , then shifts along a stroke plane inclined at an angle β. The translational and rotational motions of the hovering wing are described as follows

A(t) = A 0 2 cos( 2πt T ) + 1 , (84)
θ(t) = θ 0 1 -sin( 2πt T + φ 0 ) , (85)
where A 0 is the translational amplitude, 2θ 0 the rotational amplitude, T the flapping period and φ 0 the phase difference. The chord length c and the maximum velocity U max = πA 0 /T along the flapping path are used as the length and the velocity scales, respectively. The Reynolds number is defined as Re = U max c/ν. We employ the same parameters as used in [START_REF] Wang | Two dimensional mechanism for insect hovering[END_REF][START_REF] Xu | An immersed interface method for simulating the interaction of a fluid with moving boundaries[END_REF][START_REF] Yang | A simple and efficient direct forcing immersed boundary framework for fluid-structure interactions[END_REF] As suggested by Yang and Stern [START_REF] Yang | A simple and efficient direct forcing immersed boundary framework for fluid-structure interactions[END_REF], this simulation is performed on a large square domain of [-24c, 24c] × [-24c, 24c] to obtain a better periodicity for the results. A uniform mesh of 2400 × 2400 is employed to cover the computational domain and the mesh spacing around the wing is 0.02c, which is slightly finer than the grid resolution used in [START_REF] Xu | An immersed interface method for simulating the interaction of a fluid with moving boundaries[END_REF][START_REF] Yang | A simple and efficient direct forcing immersed boundary framework for fluid-structure interactions[END_REF]. A larger time step is selected in the present study (∆t = 0.01) based on the CFL number (CFL max = 0.72), while a much smaller time step ∆t = 0.001 is used in the immersed interface method (IIM) of Xu and Wang [START_REF] Xu | An immersed interface method for simulating the interaction of a fluid with moving boundaries[END_REF] to reduce the body shape distortion.

: c = 1, e = 4, A 0 = 2.5c, θ 0 = π/4, T = πA 0 /c, β = π/3, φ 0 = 0, Re = 157.
Figure 31 shows the vorticity fields near the flapping wing in one flapping period at four different positions, which are very similar to those given in [START_REF] Wang | Two dimensional mechanism for insect hovering[END_REF][START_REF] Xu | An immersed interface method for simulating the interaction of a fluid with moving boundaries[END_REF][START_REF] Yang | A simple and efficient direct forcing immersed boundary framework for fluid-structure interactions[END_REF]. A pair of leading and trailing edge vortices of opposite rotation is formed into a dipole. The dipole moves downward, generating the The time history of the drag and lift coefficients are plotted in Figure 32 and compared to the results of [START_REF] Wang | Two dimensional mechanism for insect hovering[END_REF][START_REF] Xu | An immersed interface method for simulating the interaction of a fluid with moving boundaries[END_REF][START_REF] Yang | A simple and efficient direct forcing immersed boundary framework for fluid-structure interactions[END_REF]. Good agreements have been found. Note that in order to maintain the shape of the rigid body in the immersed interface method of [START_REF] Xu | An immersed interface method for simulating the interaction of a fluid with moving boundaries[END_REF], a feedback control technique is employed and the time step is kept small to reduce the shape distortion. The present immersed boundary method is found to be much more satisfactory, since no additional springs for feedback control are needed and the no-slip boundary condition is exactly imposed at the interface.

A grid convergence study is also conducted to assess the accuracy of current MIBM in this case. A domain size of [-4D, 4D] × [-4D, 4D] is chosen and the grid spacing varies sequentially. The numerical solution after one flapping period is used for the analysis. A fine time step of 10 -4 is selected in order to ensure the analysis is not influenced by the temporal discretization error. Figure 33 shows the error of the horizontal velocity in L 2 norm as a function of the grid spacing. A convergence rate of around 1.29 is observed.

Flow past an impulsively started cylinder

As our last example we present results of a suddenly accelerated circular cylinder in a quiescent fluid at different Reynolds numbers Re = U 0 D/ν ranging from 40 to 3000, with U 0 being the cylinder moving velocity. Initially we place the cylinder with unit diameter (D = 1) at the origin and suddenly set it into motion to the left at a constant velocity U 0 = -1, as illustrated in Figure 34.

We first consider the Reynolds number Re = 40 and compare our results to the IBPM of Taira and Colonius [START_REF] Taira | The immersed boundary method: A projection approach[END_REF]. A uniform grid is used to cover the computational domain with no-slip boundary condition applied at all outer boundaries. The grid resolution is h = 0.01D and the time step is set to ∆t = 0.001. Two computational domains are employed to examine the effect of finite domain size on the results, namely a large domain of [-16.5D, 13.5D] × [-15D, 15D] as used by Taira and Colonius [START_REF] Taira | The immersed boundary method: A projection approach[END_REF] and a relative smaller domain [-8D, 4D] × [-5D, 5D] as used by Mimeau et al. [START_REF] Mimeau | Vortex penalization method for bluff body flows[END_REF]. The time history of the drag coefficient is plotted in Figure 36a results are in excellent agreement with the immersed boundary projection method [START_REF] Taira | The immersed boundary method: A projection approach[END_REF] on the large computational domain. When the computational domain is reduced the resulting drag coefficient is increased, which has also been observed in previous test cases. The snapshots of the vorticity field are shown in Figure 35a. Good agreements have been found compared to IBPM of Taira and Colonius [START_REF] Taira | The immersed boundary method: A projection approach[END_REF]. At this regime, a grid convergence study has been performed on a domain [-2D, 2D] × [-2D, 2D]. The time step is set to ∆t = 0.0001 and the grid spacing changes sequentially. The numerical errors are computed at t = 0.5 based on a very fine grid. Figure 37 shows the variation of the L 2 norm error for the horizontal velocity as a function of the grid spacing. A little better than first order spatial accuracy is observed.

Next we increase the Reynolds number to Re = 550 and compare our results to the vortex methods of Koumoutsakos and Leonard [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF] and Mimeau et al. [START_REF] Mimeau | Vortex penalization method for bluff body flows[END_REF]. In this case, the computational domain [-8D, 4D] × [-5D, 5D] is used and the mesh resolution is set to h = 0.005D as suggested by Mimeau et al. [START_REF] Mimeau | Vortex penalization method for bluff body flows[END_REF]. The time step ∆t = 0.001 is used. The time evolution of drag coefficient is displayed in Figure 36b. The current method has difficulties in drag prediction at early times of impulsive motion, which is also encountered by the immersed boundary projection method of Taira and Colonius [START_REF] Taira | The immersed boundary method: A projection approach[END_REF] and the vortex penalization method of Mimeau et al. [START_REF] Mimeau | Vortex penalization method for bluff body flows[END_REF]. At later stage, our results are comparable to those using vortex method. The corresponding vorticity fields are shown in Figure 35b, which compare well with the simulation results of [START_REF] Mimeau | Vortex penalization method for bluff body flows[END_REF][START_REF] Mittal | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF][START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Ploumhans | Vortex methods for high-resolution simulations of viscous flow past bluff bodes of general geometry[END_REF].

At Re = 1000, the grid is further refined to h = 0.0025D in order to solve the very thin boundary layer, while the computational domain [-8D, 4D] × [-5D, 5D] is kept unchanged. The time step is reduced to ∆t = 0.0005. As mentioned by Mimeau et al. [START_REF] Mimeau | Vortex penalization method for bluff body flows[END_REF], the two-dimensional simulation performed here is valid since only the impulsive start of the flow is considered before the onset of three-dimensional instabilities. Figure 36c and Figure 35c show the drag time evolution and the snapshots of vortex structures at different stages, respectively. We notice that the predicted drag coefficient with present method is slightly higher than that with vortex methods [START_REF] Mimeau | Vortex penalization method for bluff body flows[END_REF][START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF]. This can be attributed to the finite domain size used in the present study.

Finally we increase the Reynolds number to Re = 3000. At this Reynolds number, the simulation is quite challenging as it requires a very fine grid to capture the boundary layer. We reduce the grid size to h = 0.00125D and adjust the time step respectively to ∆t = 0.0002. Due to memory limits, we select a much smaller computational domain

CONCLUSIONS

We presented a new implicit but very efficient formulation of immersed boundary method for simulating incompressible viscous flow past complex stationary or moving boundaries. The current method treats the boundary force and the pressure as Lagrange multipliers for satisfying the no-slip and the divergence-free constraints. The fractional step method is applied to decouple the pressure as well as the boundary force from the fluid velocity field, and the two Lagrange multipliers are solved separately within their own systems. The main advantages of current approach are the accurate imposition of the no-slip condition and the efficiency in computation. The system matrices are well conditioned and generic solvers can be used directly. Especially for moving boundaries, only the boundary force coefficient matrix is updated while the coefficient matrices of velocity and pressure remain unchanged. Even though we have only dealt with rigid boundary in this article, deformable body with its motion known a priori can also be handled. A variety of distinct two dimensional flows are simulated and the results are in excellent agreement with available data sets in the literature, demonstrating the fidelity of the proposed method.

  G. CAI ET AL.

Figure 1 .

 1 Figure 1. Staggered mesh arrangement for the pressure and the velocity.

6 S

 6 .-G. CAI ET AL.

  not smooth in the support domain. Peskin[START_REF] Peskin | The immersed boundary method[END_REF] constructed a smoothed 4-pointwidth function as follows

13 ( 5 )

 135 Correct the fluid velocity with the boundary force to account for the immersed objects

  G. CAI ET AL.

F

  G. CAI ET AL.

Figure 6 .

 6 Figure 6. Global structure of the moving immersed boundary method.

  G. CAI ET AL.

Figure 7 .

 7 Figure 7. Contour of the scalar field after the boundary forcing: (a) The explicit direct forcing IBM of Uhlmann [10]; (b) The improved explicit direct forcing IBM of Kempe and Fröhlich [14]; (c) The multidirect forcing IBM of Luo et al. [15] and Breugem [16]; (d) Present MIBM.

Figure 9 .

 9 Figure 9. Temporal (a) and spatial (b) convergence analysis of current fluid solver and moving immersed boundary method for the decaying vortices problem.

Figure 10 .

 10 Figure 10. Vorticity contours and streamlines of the lid-driven cavity flow with a cylinder at Re = 1000, where the vorticity contour value is varied from -3 (blue) to 3 (red) with an increment of 0.4. Results of present MIBM are listed on the left; Results of the body-conforming mesh method are on the right.

Figure 11 .

 11 Figure 11. Comparison of velocity profiles of the lid-driven cavity flow with a cylinder at Re = 1000: (a) Distribution of velocity component u along x = 0.5; (b) Distribution of velocity component v along y = 0.5.Solid lines represent current method and dashed lines are the traditional body-conforming mesh method.

Figure 12 .

 12 Figure 12. L 2 error norm of the horizontal velocity component (u) as a function of grid spacing for the lid-driven cavity flow with an embedded cylinder.

  in present work can no longer maintain the second order accuracy. Beyer and LeVeque

Figure 13 .= 40 Figure 16 .

 134016 Figure 13. Sketch of the flow over a stationary circular cylinder.

Figure 17 .

 17 Figure 17. The wall pressure coefficient Cp and the wall vorticity Wz for flow over a stationary cylinder at Re = 40. --, results of boundary-fitted grid of Braza et al. [29]; , present h = 0.04D; △, present h = 0.029D; +, present h = 0.02D.

02D) 1 .Figure 21 .

 121 Figure 21. The wall pressure coefficient Cp and the wall vorticity Wz for flow over a stationary cylinder at Re = 100. Time-averaged values are used. --, results of boundary-fitted grid of Braza et al. [29]; , present h = 0.04D; △, present h = 0.029D; +, present h = 0.02D.

Figure 22 .

 22 Figure 22. Distribution of the horizontal velocity error on the 90 × 90 grid for the flow over a stationary circular cylinder.

Figure 23 .

 23 Figure 23. L 2 error norm of horizontal velocity u versus the computational grid size for the flow over a stationary circular cylinder.

Figure 24 .Figure 28 .

 2428 Figure 24. Time evolution of drag and lift coefficients for the flow over a stationary cylinder at Re = 1000.

Figure 29 .

 29 Figure 29. L 2 norm of the horizontal velocity component (u) versus grid spacing for the oscillating cylinder problem.

Figure 29 shows

 29 Figure29shows the results of convergence study on a domain of [-2D, 2D] × [-2D, 2D]. A time step of 10 -4 is selected and the calculation is performed for 2000 time steps. A slightly better than first order accuracy is found in this case.

Figure 30 .

 30 Figure 30. Configuration for flow over a flapping wing.

Figure 33 .

 33 Figure 33. L 2 norm error of the horizontal velocity component (u) as a function of grid spacing for the flapping wing problem.

Figure 34 .

 34 Figure 34. Sketch of the flow past an impulsively started cylinder.

  [-4D, 2D] × [-3D, 3D]. The temporal G. CAI ET AL.

Figure 35 .Figure 36 .

 3536 Figure 35. Computed vorticity contours for a suddenly started cylinder at different stages in the start-up process. Contour levels are set from -3 to 3 in increments of 0.4.

Figure 37 .

 37 Figure 37. L 2 norm error of velocity (u) for the impulsively started cylinder problem.

ACKNOWLEDGEMENT

  The first author greatly acknowledges the financial support of the China Scholarship Council. The calculations were performed on the platform PILCAM2 at the Université de Technologie de Compiègne and the HPC at the Université de Strasbourg.

  . The resulting homogeneous Neumann boundary condition of φ n+1 implies

	that	∂p n+1
	P
	e e
		r
		R
		e v
		i e w
		O n l
		y

Table I

 I

	Methods of parallelization Cores		CG	CG+MG
				Time (s) Speed-up Time (s) Speed-up
	CPU	1		63.00	1.0	25.25	1.0
		2		30.68	2.05	10.59	2.38
		4		15.26	4.13	4.56	5.23
		8		7.79	8.09	2.13	11.84
		16		4.20	15.00	1.10	22.92
		20		8.88	7.09	1.11	22.71
	GPU	240	10.36	6.08	0.63	40.08
	None	0	10.36		1.0	10.36	1.0
	AINV	1.26	6.52		1.59	7.78	1.33
	MG	0.51	0.12		86.33	0.63	16.44

. Time consummation and speed-up of the CPU and GPU parallelization for solving the pressure Poisson equation on a 400 × 400 grid. The tolerance is set to 1 × 10 -10 . Preconditioner Construction time (s) Application time (s) Speed-up Total time (s) Speed-up

Table II

 II

. Comparison of different preconditioners in GPU parallelization with the CUSP library, where the CG solver is used for solving the PPE on a 400 × 400 grid. The tolerance is set to 1 × 10 -10 .

Table III .

 III Comparison of the computational time and the velocity error. The iteration number is fixed for the explicit methods of Uhlmann[START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] and Kempe and Fröhlich[START_REF] Kempe | An improved immersed boundary mehod with direct forcing for the simulation of particle laden flows[END_REF], while others are solved until convergence under a tolerance of 1 × 10 -15 .

			Process time (s)
		Interpolation Forcing	Distribution	Total	Iter.	Error
	Uhlmann [10]	2.77 × 10 -3	1.00 × 10 -6	3.23 × 10 -3	6.02 × 10 -3	1	3.01 × 10 -1
	Kempe and Fröhlich [14]	8.15 × 10 -3	1.00 × 10 -6	8.92 × 10 -3	1.71 × 10 -2	3	7.41 × 10 -2
	Luo et al. [15] and Breugem [16] 1.16 × 10 1	1.17 × 10 -3	1.31 × 10 1	3.65 × 10 1	4443 9.96 × 10 -16
	Present	4.32 × 10 -4	1.19 × 10 -4	4.41 × 10 -4	1.33 × 10 -2	60	8.29 × 10 -16
	10 0
	10 -3
	10 -6
	10 -9
	10 -12
	10 -15	100 200 300 400 500 600 700 800 900 1000

Table VI .

 VI Comparison of the drag, lift coefficients and the Strouhal number for the flow around a stationary cylinder at Re = 100, 200. The experimental results are marked with (⋆).

		C D	C ′ D	C ′ L	St
	Re = 100 Williamson [48] ⋆	-	-	-	0.164
	Uhlmann [10] Ji et al. [18] Braza et al. [29] Liu et al. [49] Mimeau et al. [50] Xu and Wang [51] Present (Ω = 30D × 30D, h = 0.04D) Present (Ω = 30D × 30D, h = 0.029D) 1.377 ±0.010 ±0.337 0.160 1.453 ±0.011 ±0.339 0.169 1.376 ±0.010 ±0.339 0.169 1.359 ±0.019 ±0.293 0.16 1.350 ±0.012 ±0.339 0.165 1.40 0.165 ±0.010 ±0.32 0.171 1.423 ±0.013 ±0.34 1.380 ±0.010 ±0.343 0.160 Present (Ω = 30D × 30D, h = 0.02D) 1.379 ±0.010 ±0.346 0.160 Present (Ω = 40D × 40D, h = 0.029D) 1.366 ±0.010 ±0.342 0.160 Present (Ω = 60D × 60D, h = 0.029D) 1.353 ±0.010 ±0.335 0.160
	Re = 200 Williamson [48] ⋆	-	-	-	0.197
	Taira and Colonius [17] Ji et al. [18] Braza et al. [29] Liu et al. [49] Mimeau et al. [50] Xu and Wang [51] Present	1.35 1.354 ±0.044 ±0.682 0.20 0.196 ±0.048 ±0.68 1.386 ±0.040 ±0.766 0.20 1.31 0.192 ±0.049 ±0.69 1.44 ±0.05 0.200 ±0.75 1.42 ±0.04 0.202 ±0.66

Table VII

 VII

	30	S.-G. CAI ET AL.
	P
	e e
	r
	R
	e v
		i e w
		O n l
		y

Table VII .

 VII Effects of different discrete delta functions on the drag, lift coefficients and the Strouhal number for the flow around a stationary cylinder at Re = 100 and 200.

	C D	C ′ D	C ′ L	St
	Re = 100 φ 2 1.388 ±0.010 ±0.346 0.166 φ 3 1.377 ±0.010 ±0.339 0.166 φ 4 1.379 ±0.011 ±0.343 0.166
	Re = 200 φ 2 1.391 ±0.047 ±0.709 0.198 φ 3 1.365 ±0.044 ±0.696 0.200 φ 4 1.358 ±0.045 ±0.688 0.195