An efficient implicit direct forcing immersed boundary method for incompressible flows

A novel efficient implicit direct forcing immersed boundary method for incompressible flows with complex boundaries is presented. In the previous work [1], the calculation is performed on the Cartesian grid regardless of the immersed object, with a fictitious force evaluated on the Lagrangian points to mimic the presence of the physical boundaries. However the explicit direct forcing method [1] fails to accurately impose the non-slip boundary condition on the immersed interface. In the present work, the calculation is based on the implicit treatment of the artificial force while in an effective way of system iteration. The accuracy is also improved by solving the Navier-Stokes equation with the rotational incremental pressure- correction projection method of Guermond and Shen [2]. Numerical simulations performed with the proposed method are in good agreement with those in the literature.


Introduction
Immersed boundary method (IBM) is a type of non-body-conforming method, which significantly simplifies the flow simulations with arbitrary complex boundaries. Our work is based on the Uhlmann's explicit direct forcing immersed boundary method [1], which is known not to precisely satisfy the non-slip boundary condition on the interface. To improve this, we propose a low computer-consuming implicit scheme in conjunction with the rotational incremental pressure-correction projection method of Guermond and Shen [2] for solving the incompressible Navier-Stokes equations. Numerical simulations of a two-dimensional flow over a circular cylinder at two different Reynold numbers are performed, which are in good agreement with the previous experimental and numerical results in the literature.

Numerical method 2.1. Rotational incremental pressure-correction projection method
Consider the incompressible Navier-Stokes equations where u, p, Re are the non-dimensionalized velocity vector, pressure, and the Reynolds number, respectively. The equations are discretized with finite difference scheme on a staggered grid [3]. All the spatial derivatives are approximated by the second-order central difference. The secondorder Adams-Bashforth scheme is used for the temporal discretization of convective term and the Crank-Nicolson scheme for the diffusive term. This results in the following set of equations: where H, G, D, L are the discrete convection, gradient, divergence and diffusion operators, respectively. The above equations are solved with the rotational incremental pressure-correction projection method [2], in which an auxiliary velocity u * is estimated with the pressure value at previous time level in the viscous step, and the pseudo pressure φ is used to project the predicted velocity into the divergence-free field. The two sub-steps are performed as: By taking the divergence of Eq. (6) and applying the incompressibility constraint (4), the pressure corrector is obtained The pressure is then advanced by Du * Therefore, the preceding equations are solved in the order (5), (7), (6), (8) at each time step.

Direct forcing immersed boundary method
In the direct forcing immersed boundary method of Uhlmann [1], the force F n+1 is determined on the Lagrangian points such that the desired velocity U d is satisfied on the embedded boundary where the upper-case letters indicate the quantities evaluated at the Lagrangian locations, and U * is the estimated velocity by all the explicit terms in the momentum equation. Therefore, the non-slip boundary condition is not fully verified. As the Navier-Stokes equations are inherently implicit, we need to iterate the whole system such that F n+1 is calculated by the final divergencefree velocity U n+1 , which, however, is too cumbersome to perform.
In the present study, we propose to evaluate the force after the viscous step. Instead of iterating the whole system, we iterate the force within a small system coupled with the pressure poisson equation (PPE), which is often the most time-consuming part. As implied by Ji [4] that a considerable amount of computational time could be saved if the PPE is solved only once at each iteration within a time step no matter whether it is converged, since the intermediate pressure is not the true value. We follow this idea and iterate the small system in combination with GAMG (geometric algebraic multi-grid) pre-conditioner and BICGSTAB solver for the PPE by the library PETSc [5]. The modified algorithm is as follows: L (u n + u * ); Du * .
The three-points discrete delta function of Roma et al. [6] is employed for the regularization and interpolation.

Flow over a stationary circular cylinder
In the present work, the two-dimensional simulations of flow around a circular cylinder at two different Reynolds numbers Re = 30, 185 are performed, with the initial uniform velocity u ∞ = 1 and the cylinder diameter D = 1. The flow is simulated in a domain of [−10D, 40D] × [−15D, 15D]. Inlet boundary condition is set to uniform (u = u ∞ , v = 0), while ∂u/∂x = 0 is applied at the outlet plane. Free-slip conditions are assigned at lateral boundaries.
At Re = 30, the flow presents steady-state characteristics with a recirculating region in the wake of the cylinder. Whose dimensions are characterized by the length of the wake l, the distance a from the cylinder to the vortex center, the distance b between two symmetry vortex centers, and the angle θ of flow separation measured from x-axis. The wake properties in the present study is compared against previous experimental and numerical studies in Table 1 and proved to be in accord. The streamlines and vorticity contours for this flow are shown in Fig.  1. At higher Reynolds number Re = 185, periodic vortex shedding is reproduced, Fig. 1

Conclusions
An implicit direct forcing immersed boundary method is developed with a small coupled system solved in an efficient way. The rotational incremental pressure-correction projection method is employed to enforce the incompressibility constraint. The accuracy of the proposed method has been proven by the numerical simulations and comparisons with the previous experimental and numerical result in the literature.