A sharp interface direct-forcing immersed boundary method using the moving least square approximation

Based on the Moving Least Square (MLS) approximation, we propose a sharp interface direct-forcing immersed boundary method for incompressible ﬂuid ﬂows with ﬁxed and moving boundaries. Since the domain of deﬁnition for the interpolation is highly ﬂexible and the MLS approximation provides an accurate reconstructed approximation of the solution, the proposed method serves the precision and versatility required for a numerical framework to study the ﬂuid-structure interaction problems. To alleviate the inherent spurious numerical oscillation that occurs in the calculated forces on moving boundary embedded objects, we use a two step predictor-corrector method in which the direct forcing terms are calculated after the predictor step and imposed on the whole solid domain as well as at the immediate vicinity of the solid boundary inside the ﬂuid domain. To represent the arbitrary geometries, we adopt a signed distance function representation of the rigid body and an interpolation strategy to considerably reduce the computational cost of the re-initialization of the distance function at every time step. The potential capability of the method is demonstrated for both ﬁxed and moving boundary problems. We also solve a sedimentation of a single cylinder to demonstrate the ability of the present method in solving ﬂuid-structure interaction problems. These numerical experiments show that the proposed moving least square immersed boundary method can handle relatively complex moving problems while enjoying a versatile interpolation strategy and keeping the boundary conditions sharp with remarkable accuracy.


Introduction
Immersed boundary methods as a well-established sub-category of Cartesin mesh methods are powerful approaches to study complex fluid flows and fluid-structure interaction problems in engineering, cardiovascular and biological discipline, to name but a few, due to their ability to handle with relative ease of multi-body and deformable objects. Based on the strategy that is used to consider the presence of immersed objects in the flow, these methods are broadly categorized into two general classes: diffuse interface (DI) and sharp interface (SI) methods.
The DI method was originally introduced in 1972 by Peskin (Peskin, 1972). In his pioneering work, Peskin dealt with the solid boundaries by adding a singular forcing term to the momentum equations in the vicinity of the body and distributed them onto the Cartesian grid points by convoluving the field with discrete delta functions. Since those discrete delta functions smear out the force among the layers of grid points, the name of diffused interface stems from this smearing process at the surface of the solid boundary. This method was extended later in several ways; among the others, Yuki et al. (2007) introduced the solid volume fraction instead of the Lagrangian markers which also leads to a smooth interface.
In the SI methods, the immersed boundary condition is implemented at the vicinity of the Lagrangian solid boundary either by changing the background grid cells around the immersed object with the so-called cut-cell methods (Ye et al., 1999) or by using a direct forcing method in which the velocity field is reconstructed according to the desired boundary condition (Mohd-Yusof, 1997;Uhlmann, 2005;Kor et al., 2017). Deploying a local unstructured mesh, the cut-cell Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021)  domain. Liao et al. (2010) proposed this approach with a linear interpolation along the line normal to the surface of the immersed object and revealed the superiority of the method for spurious oscillation reduction.
Aiming at reducing the inherent spurious oscillation of the immersed boundary method, introduced above, and increasing the accuracy of the force calculation in direct forcing immersed boundary methods, we propose a novel approach which retains the simplicity of standard predictor-corrector direct forcing method while eliminates the need for an ad-hoc directional interpolations, by utilizing the desired feature that comes by the mesh-free interpolation. Calculation of the required forces representing sharp solid-fluid interfaces, to the best of our knowledge, has not been proposed yet with a moving least square interpolation (approximation) methodology in the framework of a predictor-corrector direct forcing method, and our numerical experiments reveals its merits for the fluid-structure interaction problems in particular for moving boundary problems. Specifically, we point out that the present method not only enjoys the moving least square approximation for data transformation between Eulerian and Lagrangian nodes (in contrast to Vanella and Balaras's (2009) approach, which results in a DI) but also it improves the quality of the solution comparing to the original direct forcing by employing a relatively simple technique to reduce the spurious oscillations. This technique closely resembles the seamless immersed boundary methods (Tajiri et al., 2014) though with a more accurate predictor-corrector strategy and interpolation which leads to a remarkable less spurious oscillation solution for moving boundary cases. This paper is organized as follows. In section 2, the governing equations for incompressible fluid flows as well as the solid body are introduced. The cell-centered finite difference discretization method used in the present study is also briefly discussed. Further, the details of the moving least square approximation is provided. We also describe the key aspect of the present two-step predictor-corrector direct forcing immersed boundary method, including the interface reconstruction and the geometrical considerations for the interpolation and its domain. In Section 3, we demonstrate the superiority of the present moving least square direct forcing formulation by applying it to the following benchmark problems: the flow past a circular cylinder at different Reynolds numbers, the flow over an oscillating cylinder in fluid at rest, and fluid-solid interaction between a falling cylinder and the surrounding fluid. Section 4 summarizes the research finding and discuss the future research direction.

Methodology of the immersed boundary technique
In this section, the general mathematical formulation of the fluid structure interaction problem is introduced. Then, the employed numerical method is discussed, followed by the details on the construction of the immersed boundary forces needed for representing the solid objects.

Governing equations
We consider the rigid body motion within an incompressible fluid, as schematically depicted in Fig. 1 with the parameters involved in this fluid-structure interaction problem. In Fig. 1, the whole domain is indicated by Ω while the fluid and the rigid body occupy the sub-domains Ω f and Ω s , respectively. The whole system is assumed to be under the gravitational acceleration g and the interaction occurs at the common boundary between the solid object and the fluid ∂Ω i = Ω f ∩ Ω s . In addition, ω s stands for the angular velocity, u s for the translational velocity and ∂Ω f for the outer boundary of the system. Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021) The fluid motion is governed by the incompressible continuity and Navier-Stokes equations, where u is the fluid velocity vector. The stress tensor σ f is given by where p is the pressure, ρ f is the fluid density, and ν is the kinematic viscosity. To make sure that the problem is wellposed, proper initial and boundary conditions are defined for each test case. In Eq. (1), the term f is added to account for the presence of the solid within the domain. Therefore, the fluid is simply simulated in the Ω = Ω f (t) ∪ Ω s (t). (See Fig.  1).
On the other hand, the rigid body motion for a two-dimensional problem is governed by the Newton-Euler equations, i.e., where m s , ρ s and I s are the solid mass, density, and the moment of inertia, respectively. The force F s and torque T s acting from the fluid on the object are obtained as where r = x s − x c is the location vector of the surface point x s with respect to the center of the mass x c . For this object, n is the outward normal vector to the surface ∂Ω i . The surface forces can be recast to the body forces using the Gauss' theorem and the immersed boundary forces (i.e., f in Eq. (1)) as (De, 2018) The instantaneous location of each rigid object is obtained by integrating the ordinary differential equations for their kinetics where θ c represents the rotation angle around the center of mass. To consider the effect of the fluid structure interaction, the boundary condition at the interface, i.e., is also imposed.
2.2. Numerical schemes 2.2.1. Overview of the numerical procedure The numerical procedure used here consists of a spatial discretization on a collocated grid and a fractional step method for the velocity-pressure coupling. In this study, all variables are evaluated at the cell centers of a Cartesian mesh, and the spatial derivatives are evaluated using the second-order central finite difference method. However, to cope with the checkerboard problem in velocity-pressure coupling, the contravariant component of velocity normal to each cell surface is also added to the list of variables which is updated during the time integration. In order to evolve Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021) the system of equations (1) in time, they are discretized using a fractional step (projection) method being implemented with a combination of the Adams-Bashforth and Crank-Nicolson schemes for the advection and diffusion terms of the momentum equation, respectively. Several variants of the fractional step methods are available, such as Kim and Moin (1985), Bell et al. (1989), and Guermond et al. (2006) with the references there in, among the others. Here we employ the formulation of Kim and Moin (1985) due to its straightforward boundary condition enforcement.
The discrete form of the algorithm entails the following sequence: where the superscripts n − 1, n, and n + 1 show the discrete times, andũ is the first provisional velocity between time step n and n + 1, which is a non-divergence-free approximation for velocity at time step n + 1. The discrete spatial derivative operators (i.e., second-order central finite difference) are denoted by ∇ h , and U is the contravariant component of the velocity vector residing at the cell faces. The pressure normalized by the density is denoted as P = p/ρ f . In this formulation, the second provisional velocity u * is computed by where ∇ h P n cc is the pressure gradient at the cell center. Then, the provisional contravariant velocity are obtained from their cell center counterpart by an interpolation function, i.e., in which F is the interpolation function from cell center variables to the cell face. Here, the interpolation is a second-order linear one. Adding the pressure gradient on the cell face to the interpolated velocity leads to a well-coupled pressurevelocity field: According to Lee et al. (2019), this treatment will also suppress the error in the kinetic energy. The pressure corrector term φ is calculated by solving the Poisson equation under the homogeneous Neumann boundary condition on all boundaries. Finally, the velocities are updated as where ∇ h φ cc and ∇ h φ f c denote the gradients of φ at the cell center and the cell face, respectively, and the pressure is updated as Having updated the fluid velocity and pressure field, the kinetics equation for the rigid body can be solved. Therefore, Eq. (3) is discretized as where F n+1 s and T n+1 s are the exerted force and torque on the solid body at time step n + 1. Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021)

Forcing strategy
The forces applied in the momentum equation (8) to represent the effect of the rigid body must be given before the procedure discussed in Sec. 2.2.1. To this end, the required force per unit mass is simply the Newtonian acceleration of all nodes near the solid body (Balaras, 2004;Uhlmann, 2005), that is where theû is an estimation of the velocity at the time step n+1 in the neighborhood of the rigid body without considering the effect of the solid boundary, and u F is the velocity at the same location while considering the presence of the body. The second term H n+1/2 represents the convective, viscous and pressure terms in momentum equations calculated in a time level between n and n + 1, It should be stressed that here the immersed boundary force is directly calculated on Cartesian mesh. Namely, u F is set as the velocity of the rigid body for all nodes inside the solid domain Ω s , while u F is interpolated from the velocity field in the fluid domain and the velocity of three close nodes on the solid boundary for the nodes just at the vicinity of the boundary ∂Ω f . This interpolated velocity along with H is used to calculate the required force. Thus, contrarily to other studies based on moving least square reconstruction (Vanella and Balaras, 2009), the velocity is not interpolated onto the surface. Furthermore, the force is not required to be interpolated from surface to the surrounding nodes smearing the boundary. Note that the this force merely affect the nodes in the solid domain and nodes very close to it and is set to zero in the fluid domain Ω f . Details on how to choose the stencil for interpolation is discussed in the subsequent section. In order to obtain the estimated field,û, a second-order explicit time step method is utilized to evolve the velocity field from time step n to the estimated one for the time step n + 1: In contrast to the Eq. (8), in Eq. (18) an explicit time step method is used for this prediction step since an explicit treatment for time integration is more efficient in terms of computational costs. Further, it has been shown (Nicolaou et al., 2014), for the corrector step with the immersed boundary force, a semi-implicit time step results in more accurate result. Therefore, while Eq. (8) is descretized using a combination of Adams-Bashforth and Crank-Nicolson schemes, we adopt the Adams-Bashforth explicit scheme for all terms in Eq. (18).
Then, we need to calculate the velocity field taking into account the boundary condition, namely u F , which is distributed around the boundary. Aiming at this, from the explicitly estimated fieldû, the velocity of the Eulerian nodes at the vicinity of the solid boundary residing inside the fluid domain are approximated by a moving least square method such that the boundary condition is satisfied at the physical solid boundary. It is worth mentioning that apart from those points inside the fluid domain at the neighbor of the physical boundary, u F inside the body is also calculated. For the grid nodes inside the solid, the estimated velocity is set at the value consistent with both translational and rotational velocities of the body. Thus, the immersed boundary force is implemented on the nodes inside the body as well as on the fluid nodes close to it.
This procedure of velocity enforcement inside the body leads to improvement of accuracy in the results, particularly in reduction of the numerical oscillation of the variables (Liao et al., 2010;Tajiri et al., 2009). However, two issues should be resolved here: (1) the forcing location of the fluid nodes and (2) the specific details of calculating u F for those nodes. Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021)

Shape representation via a signed distance function
To represent the shape of the object, a signed distance function from the set of markers put on the solid-fluid interface, denoted as Γ, is described as where x j Γ ∈ Ω s indicates the coordinate of Lagrangian markers laid on the interface, Ω s , and n j v denotes the surface-normal vector on them, as depicted in Fig. 2(a).
For regular geometries with analytical expression, this function can be obtained readily from the mathematical definition. For more complex shapes, without a closed-form expression, one can produce the function by combining the atomic shape (regular geometry) (Günther et al., 2014) or from a Standard Triangle Language (STL) representation of the surface (Luo et al., 2014) via explicitly calculating the distance from each node to the surface elements. Aiming at increasing the versatility of the method, we initialize the distance function by calculating the distance of the each Eulerian node in the domain from the nearest Lagrangian marker. When the system evolves and the location of the boundary changes, the distance function on the Eulerian nodes must be updated to reflect the new position of objects. Here, without re-initialization, we update it according to the motion of the boundary using a radial basis function interpolation as we discuss later in this section. Figure 2(a) schematically shows how the Lagrangian markers (blue circles) are distributed on the surface of the rigid body. It is worth mentioning that the distance between markers is approximately kept equal to the resolution of the background mesh. In addition, the signed distance function φ is positive if the x i is outside the object and is negative when it is inside, while φ(x i ) = 0 indicates the surface itself.
The normal vector to the surface n j v on the Lagrangian marker point can be calculated from the normal vectors on the neighboring elements n k e (k = 1, · · · , m), which can easily be computed from the coordinates, as where α k denotes the weight for each element meeting at the Lagrangian markers. Those elements in two-dimensional cases are the line segments connecting the neighbor Lagrangian markers, and their weights are considered as their length. For example, in Fig. 2(a) the vector n 1 v is obtained from n 1 e and n 2 e , while their weights are x 1 Γ − x 2 Γ and x 2 Γ − x 3 Γ , respectively.

Node stencil for the reconstruction of the velocity
In this study, the immersed boundary force is applied on the forcing nodes. As exemplified by the point F in Fig. 2(b), the forcing node are Eulerian nodes belonging to the fluid domain and have at least one neighbor in the solid domain. Associated to each forcing node, a physical boundary point, which is defined as the closest point on the fluid-solid interface, is located by where x s is the location of the physical boundary point corresponding to the forcing nodes located at x f , and φ f and n f are the distance function and the unit normal vector to the surface at the forcing point, respectively. The unit normal vector to the surface from the forcing point can be obtained from the distance function as where ∇φ is computed using the second-order central finite difference method. In Fig. 2(b), for example, the point S 1 is the physical boundary point associated with this specific forcing node F and the normal vector to the surface is illustrated by n F . It worth nothing that the location of physical boundary points calculated by Eq. (21) may or may not coincide with the Lagrangian markers that were used to obtain the distance function.
Thanks to the versatile strategy in the proposed MLS method (section 2.2.6), an approximation of the velocity can be obtained even with one stencil point. However, in this work, the stencil of velocity reconstruction consist of three to eight points depending on the geometrical constraints such as proximity of the objects and and walls. It always have at Green circles represent the stencil used for interpolation onto the black point.
least three points from the physical boundary points and at most five fluid nodes. Thus, contrasting with traditional sharp immersed boundary method, the effect of the physical boundary is taken into account with more points, which yields sharper representation of the boundary and a higher accuracy in capturing the boundary curvature.
The algorithm for finding the required stencil points is as follow.
( 1 ) For each forcing node, the first stencil point is obviously its associated physical boundary point. This point is tagged with S 1 , as schematically shown in Fig. 2 ( 2 ) Having found S 1 , we pick two auxiliary nodes, A 2 and A 3 , in the domain surrounding the forcing point along the direction of the normal vector to surface. To spot them, we first define the unit normal vector constrained on the grid, e f , at the forcing point as where e i and n i denote the i-th component of e f and n f (i = 1, 2). Then, supposing that the indices for the forcing point F is (i, j), the indices of the auxiliary nodes would be A 1 (i + e 1 , j + e 2 ), A 2 (i + e 1 , j), and A 3 (i, j + e 2 ).
( 3 ) By using Eq. (21), the location of corresponding physical boundary points to A 1 , A 2 and A 3 are sought out and tagged as S 1 , S 2 and S 3 , respectively. The points S 1 and S 2 , and S 3 are always the first three stencil points for each forcing node.
( 4 ) Two other auxiliary nodes should be chosen (i.e A 4 and A 5 ). To do so, assuming the forcing node index, (i, j), the indices for those would be A 4 (i + 2e 1 , j), A 5 (i, j + 2e 2 ),. Each of these auxiliary points might be a different forcing node, or may be inside the solid phase of the other objects (when the objects are very close to each other) thus we must assess all points in the set of auxiliary points A = {A 1 , A 2 , A 3 , A 4 , A 5 } and select the acceptable ones from the A. At this stage, however, we just record the indices of all potential stencil nodes and later we decide which one is used for the calculation. Note that A 1 and A 2 become identical in the cases where the surface perfectly perpendicular to an axis (i.e., either n 1 or n 2 is smaller than the machine epsilon). In such exceptional cases, we seek for the nodes with a different algorithm so that A 1 = (i + e 1 , j + e 2 ), A 2 = (i + e 2 , j + e 1 ), A 3 = (i − e 2 , j − e 1 ), A 4 = (i + 2e 1 , j + 2e 2 ) and A 5 = (i + 3e 1 , j + 3e 2 ).
( 5 ) Form the set of all auxiliary points, A a subset of them belonging to the fluid domain are selected as the stencil points, while those are forcing points or the nodes inside another object(s) are ruled out. Therefore, in this study, the total number of stencil points for each forcing node is not fixed and may vary between three to eight.
It should be pointed out that in order to fully exploit the proposed moving least square approximation, one needs, at least six, distinct nodes in the stencil and typically we fulfill this requirement in our algorithm. Even in the case where both auxiliary points A 2 and A 3 are forcing points and excluded from the interpolation, we still have six authentic nodes satisfying the minimum requirement for fully accurate approximation.
However, in a highly anomalous situations, the number of stencil points of very few forcing nodes might be smaller than six. Even in those cases, the algorithm is stable and the overall accuracy is not compromised because in the worst case scenario our algorithm guarantees the existence of at least three distinct stencil nodes which is enough for acceptable accuracy (Mirzaie, 2015). Moreover, these cases hardly occur; even in their occurrence, for a typical object the number of forcing points with limited stencil are much smaller than those with normal stencil (who has at least six nodes in their stencil), thus their end result is almost negligible. For example, when an object collides the wall, for very few time steps just before and after the collision, few forcing nodes may have less than six stencil points while we have much more forcing nodes around the object with full stencil nodes; therefore, the effect of those small-stencil interpolation on the calculation of immersed boundary or physical force (e.g., lift and drag) is virtually imperceptible. Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021) 2.2.5. Method for computing the signed distance for moving objects For moving objects, frequent re-initialization of the signed distance function is mandatory. To alleviate this costly procedure, Liu and Hu (2014) suggested a method for finding the signed distance without re-initialization at every time step. We adopt this method with a slight change; here we provide a brief explanation of it for the sake of completeness.
In this strategy, a local Cartesian mesh covering the object is considered, as shown in Fig. 3. For this mesh we consider a local coordinate system (x, y), which initially coincides with the global coordinate system, denoted here as (X, Y). Note that the resolution of this local mesh is not necessarily the same as the global Cartesian mesh on which the governing equations are solved. Before time-stepping, the distance function is calculated on the local Cartesian grid. During the time integration, the object is moved together with the local coordinate system. Then, the value of the signed distance function for nodes on the global coordinate system can be interpolated from the surrounding nodes of the local mesh. While in their work, Liu and Hu (2014) used a bi-linear interpolation to obtain the desired value for the global mesh, we use a localized radial basis function (Biancolini, 2017), which is more flexible and accurate. As illustrated in Fig. 3(b), to find the value of distance function for an arbitrary node on the global mesh, such as the black circle node in the figure, a block of 3 × 3 grid nodes from the local system (x, y) is selected so that the block covers the desired node on the global mesh.
By using the displacement of mass center d and the rotation matrix R of the object, we can obtain the coordinate of the global grid node in the local coordinate system. The rotation matrix for two-dimensional cases is given as where θ is the rotation angle of the object around the z axis, The global coordinate is mapped to the local coordinate by where x is the mapped coordinate of the target node in the local system attached to the object. Then, the distance function φ of the desired node in the local coordinate is calculated by a radial basis interpolation as where r i indicates the Euclidean distance between the desired node and the node i on the local system selected for interpolation, β i are the weight coefficients that must be calculated, and P is a linear polynomial which stabilizes the system (Biancolini, 2017) and for two-dimensional problems reads P(x) = p 0 + p 1 x + p 2 y. Further, ψ is the basis function. In the present study we employed an identity function for those basis functions (i.e. ψ( r i ) = r i ); therefore in two-dimensional cases we have The weight coefficients in Eq. (26) can be obtained by solving the following system: where M is the number of grid points used for this interpolation (i.e., the green points in Fig. 3(b); in this study, M = 9.
Having determined the interface, the forcing nodes and their corresponding interpolation stencil, the moving least square method for this study is introduced in the next subsection.
Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021) 2.2.6. Moving least square approximation The moving least square method is an approximation method that reconstructs a continuous function around the point on which we need the reconstructed value, via the calculation of the weighted least squares (WLS) measure biased towards the region of interest. The idea is to take an WLS formulation for an arbitrary point, move it over the domain of computation, and compute the WLS for every point individually during this move. Here we present the classical moving least squares for which the notation is adopted from Liu (2009). Given n nodes located at x j ∈ R d ( j = 1, ..., n), the WLS approximation for a continuous function u(x) is written as limited series, i.e., where p i (x) are basis functions with m terms. In the present study, we adopt the monomial basis functions, which are a function of the spatial coordinate. The bases and the coefficient matrices can be written as a(x) T = a 0 , a x , a y , a x 2 , a xy , a y 2 .
The coefficients a(x) in Eq. (29) can be obtained by minimizing the error functional, where u j = u(x j ) and u h (x j ) is its discrete approximation that we want to obtain. The error between the desired approximation and the given function is weighted using a positive weight function w, which is a function of Euclidean distance between the evaluation point and the location of the given nodes, as detailed later. The error functional can be written in the matrix form, where p 1 (x 1 ) p 2 (x 1 ) · · · p m (x 1 ) p 1 (x 2 ) p 2 (x 2 ) · · · p m (x 2 ) . . . . . . . . . . . .
and the vector of the given values is Equation (31) can be minimized by setting ∇J = 0, i.e., If the so-called momentum matrix, M = P T W P, is not singular, the value of the coefficient is obtained with Therefore, the approximated function is computed as The weight function has an important role in the formulation of the moving least square method. These weights ensure that the approximated function is globally continuous by ensuring that nodes enter and leave the compact domain smoothly, and thus the shape function compatibility is satisfied. In the literature, a wide variety of weight functions have been proposed (Levin, 2004;Li and Liu, 2002;Biancolini, 2017). Here we use the following weight cubic spline function (Biancolini, 2017) due to its compactness and straightforward formulation for two and three dimensional problems: wherer is the distance from point x made dimensionless with respect the separation distance, that is,r = x − x i /r w and r w is the maximum distance between the forcing node and the nodes in its interpolation stencil. It should be emphasized that according to Eq. (38) the shape function and approximation would be available only when the momentum matrix is non-singular. The singularity of this matrix depends on the used basis. For the linear basis (not used in this study), the sufficient condition is that there must be at least three non-co-linear nodes in the support domain for two-dimensional problems. However, for quadratic basis the condition will be more complicated, and at least six nodes are needed in the support domain. Even with more nodes in the support domain, the invertibility of the matrix may not be guaranteed. To circumvent this issue one can enlarge the supporting domain (Li and Liu, 2002), use a modified version of the linear basis (Joldes et al., 2015) or matrix triangularization (Li and Liu, 2002). Although one might solve the singularity problem by an enlarged support domain or modified basis, the system might still be ill-posed specifically when the nodes are very close to each other, which usually happens in the field of immersed boundary method. Therefore, further tricks (Wang et al., 2018) must be taken into account for these cases if a very large supporting domain is not favorable or possible. In present study, we have used a singular-value decomposition (SVD)-based pseudo-inverse method, where both singular and ill-conditioned problems are bypassed by truncating the small eigenvalue values (smaller than 10 −10 ) akin to triangularizatoin techniques. Applying the other methods would be possible in this framework and could be the subject of the future studies.

Numerical Results
We numerically carry out a benchmark test, namely the flow inside two concentric cylinders to reveal the spatial and temporal convergence rate of the proposed method. Also, flows around a circular cylinder are considered to validate the present method in steady-state and transient flows. Then, two moving boundary examples are solved. The flow around a predefined oscillation cylinder in a fluid at rest is an illustrative test to study the capability of the method for dealing with the inherent spurious oscillation in immersed boundary problems. Sedimentation of a cylinder in fluid is the last test case we consider to demonstrate the capability for fluid-structure interactions.

Flow inside two concentric cylinders
Here we simulate the flow between two concentric hollow cylinders as illustrated in Fig. 4(a). This test case is one of the rare cases bearing the exact solution (Taira and Colonius, 2007). The inner cylinder has radius r 1 = 0.5 whereas the Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021) outer is r 2 = 1. The angular velocity of the rotating cylinder, i.e. the inner one, is ω s and is changed smoothly from zero so that we can avoid the sudden initial move which might pose some adverse effects on the time accuracy of the solution, i.e., ω s = u θ (r 1 ) where the u θ is the tangential velocity component and t represents the time.
The outer cylinder is kept fixed and the initial condition is deemed a quiescent fluid. A solid wall boundary condition is imposed on the outer boundaries of the computational domain whose size is [0, 2.1] × [0, 2.1]. We compute the azimuthal velocity along a radial line u θi (i = 1, · · · , N in r 1 ≤ r ≤ r 2 ) and report the convergence rate of its L 2 and L ∞ error norms, i.e., where the analytical solution U θ (r) is given by Figure 4(b) shows the spatial convergence rate which is slightly smaller than two for both L 2 and L ∞ . To obtain this rate, we fixed the time step at ∆t = 4.0 × 10 −6 and varied space resolution uniformly from ∆x = ∆y = 1.6 × 10 −2 to 0.41 × 10 −2 in five steps. According to this figure, the rate of decay for the spatial errors to be 1.5 and about 1.8 for L ∞ and the L 2 , respectively. Although this rate of convergence is better than diffusive direct forcing methods (Uhlmann, 2005;Taira and Colonius, 2007), which use discrete delta functions to transfer data, further investigation is required to explain why the L 2 -norm of the error is smaller than two. The impact of the splitting error in the fractional step method as well as time integration is demonstrated in Fig. 4(c) by calculating the temporal convergence rate for a spatial resolution of ∆x = ∆y = 0.41 × 10 −2 , which confirms that the temporal convergence rate is first order as designed.

Flow over a fixed circular cylinder
Next, a flow over a fixed circular cylinder is considered as a canonical test case, since a great amount of numerical and experimental studies are available for comparison. The Reynolds number Re = u ∞ D/ν based on the inflow velocity u ∞ , the cylinder diameter D, and the kinematic viscosity ν determines the flow characters. The simulation is performed in a rectangular domain with a fluid flow from left to right (see, Fig. 5). The uniform flow velocity u ∞ = 1 in the streamwise (x) direction is imposed at the inlet, the slip boundary condition is imposed at the bottom and top boundaries, and the following convective outlet boundary condition is adopted at the outlet: For the pressure, a Neumann boundary condition, ∂p/∂n = 0, is imposed on all boundaries. The cylinder is placed 20D downstream of the inlet. While a uniform coarse mesh, ∆x = ∆y = 0.04, is employed for the whole domain, the grid is refined in a rectangular domain with [−10D, 10D] × [−10D, 10D] around the cylinder so that ∆x = ∆y = h as shown in Table 1. Note that the coarse and refined meshes are smoothly connected at their interfaces by placing computational cells of intermediate sizes.
For the sake of comparison, the drag and lift coefficients are defined as where F L and F D are the exerted lift and drag force form the fluid, respectively, which can be calculated using the immersed boundary forces in Eq. (4) both on the solid domain and forcing nodes. Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021)

Steady case at Re = 40
The flow at Reynolds numbers Re = 40 is a steady flow having a re-circulation region in the wake. This region is characterized by the length of the circulation l, the streamwise distance a from the nearest point on the body surface to the vortex center, the transverse distance b between the pair of vortices, and the angle of separation θ, as depicted in Fig. 5. We perform the simulation with various domain sizes as well as different mesh resolutions to check the grid convergence and ensure that the boundary confinement effect does not impact the solution. The time step is selected small enough (∆t = 0.001) to satisfy the CFL condition (i.e., the Courant number, C < 1). Our numerical experiments revealed that the present method remains stable even with the Courant number close to one. Figure 6 shows the pressure and vorticity contours and the streamlines, which are observed to be in close agreement with those reported in the literature (Uhlmann, 2005;Taira and Colonius 2007). Table 1 shows the comparison of the wake dimension (i.e., l, a, b, and θ defined in Fig. 5) and the drag coefficient C D against those reported in the experimental and other numerical studies. The result shows a good agreement with those reference values. Also, one could conclude that the narrow domain leads to a slightly larger drag coefficient, which was also observed in the previous studies (Uhlmann, 2005;Taira and Colonius 2007).

Periodic vortex shedding cases at Re = 100 and 200
The flow over a fixed cylinder becomes unsteady and a periodic vortex shedding emerges when the Reynolds number is increased further, as can be observed in Fig. 7. In Fig. 8, the time evolution of the drag and lift for those Reynolds number are plotted. Note that the oscillating frequency of the lift is indeed the vortex shedding frequency, f s , and the Table 1 Comparison of the steady-state wake dimension and the force coefficient for flow over a circular cylinder at Re = 40. The results with * are experimental ones. frequency of the drag is twice that of the lift coefficient (Balaras, 2004). Table 2 summarizes the mean drag coefficient C D , the amplitude of drag coefficient C D and lift coefficient C L as well as the the Strouhal number, St = f s D/u ∞ . For the sake of a better comparison, the result of the experimental and numerical studies have been provided. In the numerical section, both body-fitted and immersed boundary methods are considered. A good agreement has been obtained against those well-established results. With the reference values from a body-fitted method published by Liu et al. (1998), the mean drag coefficient shows approximately merely 3% over-prediction for the smaller domain while other IBM methods such as Uhlmann (2005) overpredicts it around 11%. From Table 2, it can be seen that the accuracy will be improved by enlarging the computational domain. Further, since our method is not iterative unlike the methods of Breugem (2012) and Ji et al. (2012), the present method can be much more efficient in terms of computational costs.

In-line oscillating cylinder in fluid at rest
A predefined oscillatory motion of a cylinder in a fluid at rest, that has been frequently reported in the literature, was studied with both experimental (Dütsch et al., 1998) and numerical methods (Kor, 2017;Kor et al., 2018;Liao et al., 2010). In particular, this test case can reveal spurious oscillations in the force imposed by the immersed boundary method and is utilized to investigate the effect of treatments for those oscillations (Seo and Mittal, 2013).
The key dimensionless parameters in this flow are the Reynolds number, Re = U max D/ν and the Keulegan-Carpenter number, KC = U max /( f D), where U max and f are maximum velocity of the cylinder and the characteristic frequency of it. In the present work, we utilized the set up employed in experimental and numerical study with a non-inertial reference Table 2 Comparison of the Struhal number and force coefficients for the flow over the circular cylinder at Re = 100 and Re = 200. The results with * has been obtained from experimental study. Uhlmann (2005) 1.453 ±0.011 ±0.339 0.169 Ji et al. (2012) 1.376 ±0.010 ±0.339 0.169 Braza et al. (1986) 1.359 ±0.019 ±0.293 0.16 Liu et al. (1998) 1.350 ±0.012 ±0.339 0.165 Mimeau et al. (2015) 1.40 ±0.010 ±0.32 0.165 Xu and Wang (2006) 1  196 Ji et al. (2012) 1.354 ±0.044 ±0.682 0.20 Braza et al. (1986) 1.386 ±0.040 ±0.766 0.20 Liu et al. (1998) 1.31 ±0.049 ±0.69 0.192 Mimeau et al. (2015) 1.44 ±0.05 ±0.75 0.200 Xu and Wang (2006) 1 reported by Dütsch et al. (1998). The motion of the cylinder is predefined by its center location X C , as The Reynolds number and Keulegan-Carpenter number are set at Re = 100 and KC = 5, respectively. The computational domain has a size of 40D 40D in x and y direction, and the cylinder is initially located at the center. A uniform grid with resolution of ∆x = ∆y = 0.025 is used to discretize the domain. An appropriate time step is chosen to keep the Courant number about C 0.55 (Liao et al., 2010). Started from the static condition, the simulation continues until periodic vortex shedding is fully developed. Figure 9 shows the pressure and vorticity contours at four different phase angle (ϕ = 2π f t = 0 • , 96 • , 192 • , 288 • ) after vortex shedding begins. The vortices contour are drawn from −3 to 3 with an increment of 0.4 displaying the same structure as Dütsch et al. (1998). As the cylinder moves, upper and lower boundary layers are developed at the leading face of the cylinder and separated while the cylinder moves through them. When cylinder goes in the reverse direction, it splits the already stirred flow (i.e., vortices remaining in the field). These are consistent with the experimental observation (Dütsch et al., 1998).
In Fig. 10, the profile of velocity components are depicted in three different angle phase (ϕ = 2π f t = 180 • , 210 • , 330 • ) at four different streamwise location (x = −0.6D, 0D, 0.6D, 1.2D). The experimental results from Dütsch et al. (1998) are also plotted for comparison. The present results agree well with those from Dütsch et al. (1998) in the fluid domain. It is worth mentioning that according to this figure, the forcing treatment inside the solid phase does not influence the fluid phase and just provides a smooth transition between the fluid and the solid domains. Figure 11 shows the evolution of the drag coefficient (C D = F x /[(1/2)ρU 2 max D]) versus dimensionless time (with the period of the oscillation) at the grid resolution of ∆x = ∆y = 0.02. In this figure, the present method is compared with reconstruction direct forcing methods, and the body-fitted numerical result of the Guilmineau and Queutey (2002). The proposed method shows a much smoother drag coefficient evolution with a better agreement with the body-fitted method than the original direct forcing methods. In fact, this is the great merit of the present method, which motivated us to develop it.

A freely falling circular cylinder in the fluid at rest
Finally, the dynamics of free falling solid circular cylinder is simulated to demonstrate the ability of the present method to solve two-way fluid-structure interaction without pre-defined velocity. In this regard, the falling of a twodimensional cylinder in a rectangular tank filled with a Newtonian fluid is considered. The configuration of this test case  (Dütsch et al., 1998) is represented by red delta, green gradient, blue square, and orange circle, respectively.
is same as those of Glowinski et al. (2001). That means, the computational domain is a rectangle, (2 cm × 6 cm) , and the diameter of the cylinder is D = 0.25 cm while the initial location for the center of the cylinder center is (x, y) = (1 cm, 4 cm). The density of the fluid is ρ f = 1 g/cm 3 and that of the cylinder is ρ s = 1.25 g/cm 3 . We also take the viscosity as ν = 0.1 cm 2 /s. The mesh sizes are 1/36 and 1/64 of the cylinder diameter, that is ∆x = ∆y = (1/144) cm and ∆x = ∆y = (1/256) cm . For all of these cases, 300 Lagrangian markers are considered on the interface; thus, the geometric effects are similar for all cases. By solving equations (8)- (12) and (16), the location of the cylinder is updated and the corresponding immersed Fig. 11 The comparison of the time evolution of the drag force for a cylinder oscillating in a fluid at rest. The solid black lines obtained with the present method, the red line is obtained from the reconstruction method . The orange circular symbols indicate the body-fitted result of Guilmineau and Queutey (2002). boundary forces at each time step is calculated. We also considered a modeled collision force between cylinder and wall to prevent it from interpenetration to the wall. Aiming at keeping the physical distance between objects, Glowinski et al. (2001) proposed a short-range repulsive force model. Singh et al. (2003) modified that model by allowing a little overlap between objects. Combining these two models, Wan and Turek (2007) suggested their own repulsive model, which is utilized in the present study, i.e., where R i is the radius of cylinder, X i is its center, X i is the coordinate vector of the imaginary cylinder located on the physical boundary, and d i, j = |X i − X i |. The parameter ξ is the range of the repulsive force, while and represent the small positive parameters for the collision. For the results provided here we took = (∆x) 2 , = ∆x, and ξ = ∆x/10. Further, the maximum Reynolds number is defined as where u s (t) and v s (t) are the x and y components of the cylinder velocity. The comparison of the maximum Reynolds number during the settling in the fluid with previous numerical studies is given in Table 3. The maximum Reynolds numbers predicted under different grid resolution agree well with those of the other well-established methods (Glowinski et al, 2001;Wan and Turek, 2007). Table 4 shows the maximum Reynolds number obtained at a fixed resolution, ∆x = (1/256) cm, with different time step size ∆t, suggesting that the convergence rate for ∆t is approximately one.
The vorticity contour at four different time steps in Fig. 12 show that the flow field is symmetric and virtually no rotation happens. Also, the vorticity contours are in good agreement with the previous numerical results (Glowinski et al., 2001).  Figure 13 quantitatively compares the vertical and horizontal components for the center of the cylinder as well as the vertical location of it during the settlement with two different resolution. These results indicate that the numerical solution are well-converged and are in good agreement with those reported by the other researchers (Glowinski et al., 2001;Kor et al., 2018). In addition, Fig. 14 shows two other important quantities, i.e., the transnational kinetic energy E t = (1/2)M s (u(t) 2 + v(t) 2 )and the rotational kinetic energy E r = (1/2)I s ω s (t) 2 , where M s and I s are the mass and the momentum of the inertia of the cylinder, and ω s (t) stands for instantaneous angular velocity. It can be seen that the result computed on two different resolutions are essentially the same except for the rotational kinetic energy due to its very small amount (O(10 −4 )), which is attributed to the numerical errors.

Flow induced by a hovering elliptical wing
The hovering of an elliptical wing in a fluid initially at rest is considered to demonstrate the capability of the proposed algorithm for the cases with rotational movement. The wing undergoes a prescribed motion determined by the following horizontal movement of the center of mass and the rotation of the body: where X(t) is the time dependent location of the mass center in x direction and α is the counter-clockwise angle of the wing with respect to the positive x axis. For the sake of comparison, the parameters are set similar to Eldredge (2007): the initial angle is α 0 = π/4, the rotational amplitude is α max = π/4, and the translational amplitude is X max = 2.8 c, where c denotes the chord length of the wing. The period of the oscillation is set at T = 1/ f = 4 and the maximum translational velocity is 8.8 c/T . The Reynolds number is set at Re = 75 based on the maximum velocity and the chord length. Also, the y coordinate of the mass center is kept zero (Y(t) = 0), and the ratio of the chord length to short axis of the wing (c/b) is 10. No mean flow is assumed. Figure 15 shows snapshots of the vorticity field during 2 ≤ t/T ≤ 4. In comparison with other studies such as Eldredge (2007), our results indicate that the proposed algorithm can properly cope with problems with moving and rotating boundary. It should be noted that Eldredge (2007) obtained the results by a vortex particle method with the same configuration. Ghomizad, Kor and Fukagata, Journal of Fluid Science and Technology, Vol.16, No.2 (2021) Figure 16 compares the computed lift and drag coefficients with those of Eldredge (2007). The present results are obtained with ∆x = ∆y = 0.009 using the proposed algorithm, which show a good agreement with the reference despite the differences in the methods.

Conclusions
In this paper, we have proposed a more accurate and relatively simple sharp interface direct forcing immersed boundary method based on a two-step predictor-corrector projection strategy for solving the system of the governing equations. For the consistency consideration, we adopted the second order in time Kim-Moin scheme for the projection in the proposed method. Regarding the direct immersed boundary force term, it is calculated after prediction step, in which a time step of the projection method is taken without considering the immersed object, and distributed on the whole solid domain while the data is transferred between the Eulerian and Lagrangian nodes by using a direct moving least square approximation. In contrast to the traditional moving least square immersed boundary methods, this approach succeeded in keeping the boundary conditions sharp while reducing the spurious numerical oscillations.
We have performed several numerical experiments to illustrate the accuracy and capability of the proposed method, including flow induced by rotating cylinder inside a larger stationary hollow cylinder for convergence test, the flow past a stationary cylinder with different Reynolds number, the flow over an oscillating cylinder with predefined motion and finally the sedimentation of a cylinder due to its weight in a fluid at rest. In particular, our numerical results shows that the convergence rate of the velocity field, being close to the second order, is apparently better than those of the first order projection methods. It has been showed that our numerical is capable enough to be considered as an accurate and relatively simple method for the fluid structure interaction problem and has high potential for the other field such as heat transfer and two phase problems. As this method utilizes a versatile stencil which is not defined along any arbitrary line, it can be straightforwardly extended to the three dimensional problems. Also, the proposed method takes into account more than one nodes on the solid-fluid interface, which yields more accurate boundary capturing capability, particularly in three dimensional problems. Extension to three dimensional problem is one of our future research plans. Besides, the cell center strategy adopted in this work and the implicit interface representation pave the way for the future works on more efficient methods such as adaptive mesh refinement solvers equipped with immersed boundary method for fluid solid interaction problems.