An Immersed Boundary Method Based Improved Divergence-Free-Condition Compensated Coupled Framework for Solving the Flow–Particle Interactions

A flow–particle interaction solver was developed in this study. For the basic flow solver, an improved divergence-free-condition compensated coupled (IDFC2) framework was employed to predict the velocity and pressure field. In order to model the effect of solid particles, the differentially interpolated direct forcing immersed boundary (DIIB) method was incorporated with the IDFC2 framework, while the equation of motion was solved to predict the displacement, rotation and velocity of the particle. The hydrodynamic force and torque which appeared in the equations of motion were directly evaluated by fluid velocity and pressure, so as to eliminate the instability problem of the density ratio close to 1. In order to effectively evaluate the drag/lift forces acting on the particle, an interpolated kernel function was introduced. The present results will be compared with the benchmark solutions to validate the present flow–particle interaction solver.


Introduction
The offshore wind farm is a wind turbine installed in an offshore area where there is a strong wind field for a long period of time. It can avoid the noise, shading and visual obstruction caused by the wind turbine's operation in the land area. Offshore wind farms can also provide more stable wind energy and low turbulence effects. However, it requires consideration of the effects of wind turbulence on the wind turbine and its supporting structure. Recently, attention has also been paid to the offshore floating wind turbine. How to efficiently simulate the interaction between the floating structure and the water/air interface becomes an important research topic.
The coupling analysis of fluids and solids in practical engineering applications usually requires the consideration of complex geometries. These complex geometries are usually in a stationary or transient motion at high Reynolds number flows, such as offshore stationary wind turbines and marine systems. These fluid-structure problems are typically solved by the traditional body-fitted-grid. In order to generate high-quality mesh, a lot of efforts in terms of manpower are always required to refine the mesh at the high curvature region. The computational cost is then increased, especially when dealing with the moving boundary problems. In addition, the quality of the grid is also very crucial for simulation analysis.

Incompressible Navier-Stokes Equations
The non-dimensional primitive-variable Navier-Stokes equations for incompressible flows can be expressed as the following momentum equations and continuity equation: where u is the velocity vector, p is the pressure field, and Re is the Reynolds number. f is external force vectors due to particle motions, and is modeled by the present immersed boundary based flow-particle interaction framework.

Equations of Motion for Solid Particle
The equations of motion for the rigid body can be written as d Ω dt = ω (6) where the linear velocity vector ( V) and particle position vector ( X) are due to linear motion, while ω and Ω are the angular velocity and orientation due to the angular motion. ρ S and ρ f are the density for the solid particle and fluid, respectively. V s is the volume of the solid particle, g is the gravity, and I S is the inertial moment of the solid particle. F f and ψ f are the hydrodynamic force and torque acting on the solid particle, respectively.
As stated in [2,3], when we evaluated the F f and ψ f by the summation of forcing terms f together with the relation of momentum balance over the corresponding fluid domain [2], it will lead to an instability issue when the density ratio close to 1, or (ρ S − ρ f ) ≈ 0. This is due to the fact that the term (ρ S − ρ f ) will appear in the denominator when solving the equations of motion. In order to avoid this issue, hydrodynamic force and torque are directly evaluated in this study: ψ f = ρ f r × (σ·n) dA (8) It is noted that by using the present approach, (ρ S − ρ f ) will be never appear in the denominator, so that there will be no instability problems even when the density ratio id close to 1. In Section 3.1, the present study will introduce methods to accurately and effectively evaluate the two above equations. Furthermore, the implicit midpoint temporal scheme was employed to solve Equations (3)- (6) in the present study.

Differentially Interpolated Direct Forcing Immersed Boundary (DIIB) Method
In this section, the DIIB method is employed to model the solid body. Here, we will briefly describe the fundamental idea of the DIIB method.
Assuming there is a solid object S immersed in the Eulerian grid; in order to model the effect of solid object, the direct forcing term f is introduced in Equation (1) to satisfy the velocity condition ( V) of the solid object. With the semi-implicit discretization, Equation (1) can be split into the following two equations: It is noted that in practice, the boundary points of the solid object do not always fall into the Eulerian grid points. The direct forcing method will then evaluate the modified velocity on some target Eulerian grid points, followed by obtaining the corresponding momentum forcing terms based on Equation (10) to model the effect of solid object with velocity condition V.
For the DIIB method, the velocity (u S ) at Eulerian grid point which is inside and close to the solid boundary, will be modified to satisfy the velocity condition u IB at the solid boundary by virtue of interpolated velocity at fluid domain u f and Taylor series expansion: Taking Figure 1 as an example, in general u F does not always lie on the grid point. The evaluation of u f is done by "advect" u f from point Q to point A [18]: In this section, the DIIB method is employed to model the solid body. Here, we will briefly describe the fundamental idea of the DIIB method.
Assuming there is a solid object immersed in the Eulerian grid; in order to model the effect of solid object, the direct forcing term is introduced in Equation (1) to satisfy the velocity condition ( ⃑ ) of the solid object. With the semi-implicit discretization, Equation (1) can be split into the following two equations: It is noted that in practice, the boundary points of the solid object do not always fall into the Eulerian grid points. The direct forcing method will then evaluate the modified velocity on some target Eulerian grid points, followed by obtaining the corresponding momentum forcing terms based on Equation (10) to model the effect of solid object with velocity condition ⃑ .
For the DIIB method, the velocity ( ) at Eulerian grid point which is inside and close to the solid boundary, will be modified to satisfy the velocity condition at the solid boundary by virtue of interpolated velocity at fluid domain and Taylor series expansion: Taking Figure 1 as an example, in general does not always lie on the grid point. The evaluation of is done by "advect" from point Q to point A [18]: In the above, ⃑ is the unit normal vector of the solid object. In practice, Equation (12) is solved by first order upwind with single artificial time step ∆ = ̅̅̅̅̅̅ /|⃑ | = 2 ̅̅̅̅ /|⃑ |. The reader can refer to [18] for details. To alleviate the spurious force oscillations (SFOs) and to obtain smooth pressure, the velocity modification method [17] was employed in this study. The idea behind this method is to modify the velocity inside the velocity to the solid object moving velocity when solving the pressure correction equation. It has been shown in [17] that this method successfully suppresses the SFO for the problems with prescribed motions and will then be implemented for dealing with the problem for which the motions are obtained by solving equations of motion. In the above, n is the unit normal vector of the solid object. In practice, Equation (12) is solved by first order upwind with single artificial time step ∆τ = APQ/ n = 2AP/ n . The reader can refer to [18] for details.
To alleviate the spurious force oscillations (SFOs) and to obtain smooth pressure, the velocity modification method [17] was employed in this study. The idea behind this method is to modify the velocity inside the velocity to the solid object moving velocity when solving the pressure correction equation. It has been shown in [17] that this method successfully suppresses the SFO for the problems with prescribed motions and will then be implemented for dealing with the problem for which the motions are obtained by solving equations of motion.

Evaluation of Drag and Lift Forces
When employing the direct forcing immersed boundary method, it is suggested to evaluate forces F by the following equation for calculating drag and lift forces: where σ is the summation of shear and pressure stress, n is unit normal vector and A is surface area. In this study, a four-point piecewise discrete delta function based interpolated kernel is employed [17]: when the solid particle is approaching the wall, the interpolated kernel will be simplified to the following one-point piecewise discrete delta function [19]:

Improved Divergence-Free-Condition Compensated Coupled (IDFC 2 ) Framework
The idea behind the IDFC 2 framework is to discretize momentum equations by virtue of cell-center and cell-face velocity [17]. With the present approach, we can not only obtain the accurate velocity vector (u, v), but also the fully coupled pressure field p. The following will first introduce the original IDFC framework, then briefly describe how to mitigate to IDFC 2 framework.

Derivation of IDFC Framework
With the 2D grid structure as shown in Figure 2, the derivation of the IDFC framework starts from the semi-discrete momentum equation. Without loss of generality, we take velocity u at the cell-center P (u p ) as an example: The idea for the IDFC method is that discretizing the momentum eq cell-face should have the same form as we discretize the momentum e cell-center:  (20), * can then be obtain can also be obtained in a similar way. In order to fulfill the divergence-free Poisson equation for pressure correction is then solved with cell-face veloc It is followed by updating the cell-face velocity by virtue of the predicte rection: The idea for the IDFC method is that discretizing the momentum equation for the cell-face u e should have the same form as we discretize the momentum equation at the cell-center: However, obtaining the velocity derivatives on the cell-face is not as easy as on the cell-center. The linear-averaged approximation is then employed to evaluate the cell-face velocity derivatives on the cell-face e [20]: 1 By using Equations (21)-(23) into Equation (20), u * e can then be obtained. u * w , v * n , v * s can also be obtained in a similar way. In order to fulfill the divergence-free condition, the Poisson equation for pressure correction is then solved with cell-face velocity: It is followed by updating the cell-face velocity by virtue of the predicted pressure correction: while the pressure is updated as p n+1 = p n + p . Finally, the cell-center velocity is updated by again employing the linear-averaged approximation for the cell-center derivatives: where: It is noted that by virtue of Equations (30)-(32), the introduced IDFC forcing terms can efficiently stabilize the calculations [20].

Derivation of IDFC 2 Framework
It is known from [17] that the IDFC method may lead to non-physical velocity oscillations due to the fact that the cell-face and cell-center velocity is not strongly coupled when evaluating the velocity derivatives. To resolve this issue, we modified Equation (23) to include the contribution of cell-face velocity: In the above, C is the newly introduced parameter which is utilized as coupling cell and face velocity. In this study, C was chosen as 1/2 to ensure the strong coupling.
The approximation equation for cell-face velocity u * e is then rewritten as The corresponding cell-center velocity is expressed as The overall algorithm for the present DIIB-IDFC 2 framework is summarized as follows. The linear and angular velocity, position and rotation of the solid objects are first solved by the equations of motion, where the hydrodynamic force and torque are obtained from the fluid velocity and pressure. To compute the summation of the hydrodynamic force and torque, the solid surface is first divided by a collection of Lagrangian points, and then interpolates the fluid velocity and pressure on each Lagrangian point. Finally, the linear and angular velocity of solid objects are then employed as the boundary condition for the solid-fluid interface. For the sake of completeness, we also plotted the algorithm in Figure 3. In this study, a finite volume-based, second order accurate dispersion-relation preserving upwinding scheme [20] was utilized to discretize the convection terms, while the central difference was employed for other derivative terms. The second order semi-implicit Gear method was employed as the temporal scheme. For more details of the numerical scheme, people can refer to [20].
Energies 2021, 14, x FOR PEER REVIEW Figure 3. Overall algorithm for the present framework.

Taylor-Couette Flow
The 2-D Taylor-Couette flow problem has been chosen as the first va lem. For the present problem, the radius for outer circular ( ) and inner ci 0.4 and 0.2. The non-slip boundary condition is employed for the outer constant rotating angular velocity = 3 is set as a boundary conditio circular. The exact solution for this kind of flow setting can be derived as For the present study, four different mesh sizes, namely Δ = Δ = Δℎ 1/40 and 1/80 with = 500 are conducted for the validation study in a main. The inner and outer cylinder are modeled with DIIB method, as shown

Taylor-Couette Flow
The 2-D Taylor-Couette flow problem has been chosen as the first validation problem. For the present problem, the radius for outer circular (R O ) and inner circular (R I ) are 0.4 and 0.2. The non-slip boundary condition is employed for the outer circular, and a constant rotating angular velocity ω TC = 3 is set as a boundary condition for the inner circular. The exact solution for this kind of flow setting can be derived as For the present study, four different mesh sizes, namely ∆x = ∆y = ∆h = 1/10, 1/20, 1/40 and 1/80 with Re = 500 are conducted for the validation study in a unit square domain. The inner and outer cylinder are modeled with DIIB method, as shown in Figure 4. For the present study, four different mesh sizes, namely Δ = 1/40 and 1/80 with = 500 are conducted for the validation stu main. The inner and outer cylinder are modeled with DIIB method, a  From Table 1 and Figure 5, this shows that the rate of convergence for u, v and p are 2.033, 2.033 and 2.044, which match very well with the proposed accuracy order. The validation of the present solver was then confirmed. From Table 1 and Figure 5, this shows that the rate of convergence for , and are 2.033, 2.033 and 2.044, which match very well with the proposed accuracy order. The validation of the present solver was then confirmed.

Lid-Driven Semi-Circular Cavity Flow
To further understand the performance of the present framework, the lid-driven semi-circular cavity flow problem was then investigated. In a 1 × 0.5 rectangular domain, there is a semi-circular cavity with a radius of 0.5 modeled by immersed boundary method, and the lid is driven with velocity = 1, as shown in Figure 6. Comparisons of cutting-line velocity at = 0.5 and = 0.25 with the benchmarking solutions of Glowinski et al. [21] and Ding et al. [22] were made with = 1000 and 100 × 50 mesh. This shows excellent agreement with Figure 7. To further show the accuracy of the present IDFC 2 framework, we also plotted the numerical results which was obtained by the fifth order upwinding scheme-based DFC framework [23]. As can be seen in Figure 7, the present results match better than the reference results. For the sake of completeness, the predicted contours of pressure and vorticity are also plotted in Figure 8 to show the applicability of the present framework.

Lid-Driven Semi-Circular Cavity Flow
To further understand the performance of the present framework, the lid-driven semicircular cavity flow problem was then investigated. In a 1 × 0.5 rectangular domain, there is a semi-circular cavity with a radius of 0.5 modeled by immersed boundary method, and the lid is driven with velocity u lid = 1, as shown in Figure 6. Comparisons of cutting-line velocity at x = 0.5 and y = 0.25 with the benchmarking solutions of Glowinski et al. [21] and Ding et al. [22] were made with Re = 1000 and 100 × 50 mesh. This shows excellent agreement with Figure 7. To further show the accuracy of the present IDFC 2 framework, we also plotted the numerical results which was obtained by the fifth order upwinding schemebased DFC framework [23]. As can be seen in Figure 7, the present results match better than the reference results. For the sake of completeness, the predicted contours of pressure and vorticity are also plotted in Figure 8 to show the applicability of the present framework.
there is a semi-circular cavity with a radius of 0.5 modeled by immersed boundary method, and the lid is driven with velocity = 1, as shown in Figure 6. Comparisons of cutting-line velocity at = 0.5 and = 0.25 with the benchmarking solutions of Glowinski et al. [21] and Ding et al. [22] were made with = 1000 and 100 × 50 mesh. This shows excellent agreement with Figure 7. To further show the accuracy of the present IDFC 2 framework, we also plotted the numerical results which was obtained by the fifth order upwinding scheme-based DFC framework [23]. As can be seen in Figure 7, the present results match better than the reference results. For the sake of completeness, the predicted contours of pressure and vorticity are also plotted in Figure 8 to show the applicability of the present framework.

Flow Past Circular Cylinders in Tandem
In this subsection, the flow past circular cylinders in tandem

Flow Past Circular Cylinders in Tandem
In this subsection, the flow past circular cylinders in tandem problem will be simulated to show that the present framework can produce accurate solutions for time transient flows. For this problem the computational domain is set as a 45 × 20 rectangular domain at the inlet with velocity ∞ = 1 from the left boundary, while the other three boundaries are set as convective boundary conditions. Inside the computational domain, there are two circular cylinders with diameter = 1 in tandem. The centroid of the left one is located at (9.5, 0), while the right one is located at (9.5 + D, 0), where D is the distance between two circular centroids. In this study, two scenarios of D = 4 and D = 5 with =

Flow Past Circular Cylinders in Tandem
In this subsection, the flow past circular cylinders in tandem problem will be simulated to show that the present framework can produce accurate solutions for time transient flows. For this problem the computational domain is set as a 45 × 20 rectangular domain at the inlet with velocity u ∞ = 1 from the left boundary, while the other three boundaries are set as convective boundary conditions. Inside the computational domain, there are two circular cylinders with diameter d = 1 in tandem. The centroid of the left one is located at (9.5, 0), while the right one is located at (9.5 + D, 0), where D is the distance between two circular centroids. In this study, two scenarios of D = 4 and D = 5 with Re = 200 and mesh size ∆x = ∆y = d/20 were investigated. The details of the present settings are plotted in Figure 9. The simulated vorticity and pressure contours at t = 400 are plotted in Figures 10 and 11 for D = 4 and D = 5, respectively. As shown in Figures 10 and 11, the present framework can correctly produce the vortex shedding behavior, while the pressure can remain smooth. The predicted results were then used to evaluate the drag coefficient (C D ), lift coefficient (C L ) and Strouhal number (St), and compare with the referenced solutions ( [24][25][26][27]). From Table 2, good agreements can be observed between predicted results and benchmarking solutions. The drag and lift coefficient are also plotted in Figures 12 and 13

Two Circular Cylinders Moving towards Each Other in Quiescent Flow
In the computational domain ( = −8~24, = −8~8), there are two circular cylinders with a unit diameter (d = 1), the lower initially located at ( , = 0,0), while the upper

Two Circular Cylinders Moving towards Each Other in Quiescent Flow
In the computational domain ( x = −8 ∼ 24, y = −8 ∼ 8), there are two circular cylinders with a unit diameter (d = 1), the lower initially located at (x l , y l = 0, 0), while the upper one located at (x u , y u = 16, 1.5), and then moving toward each other along the x direction with the prescribed motions until t = 32: Based on the above prescribed motions, two cylinders will move periodically in the x axis until t = 16, then start to move towards each other. In this study, a no-slip boundary condition is applied to all boundaries, and the mesh size and time step size are chosen as ∆x = ∆y = d/40 and ∆t = 10 −2 . From Figures 14 and 15, it shows that the present framework can produce smooth pressure and vorticity. To validate the accuracy of the predicted results, we plotted Figure 16 to compare the time-history drag (C D ) and lift (C L ) coefficients of the upper cylinder with the referenced solutions [9,28] which shows excellent agreements. For the sake of completeness, we also plotted Figure 17 to show the effect of velocity modification method. This clearly shows that with the velocity modification method, SFO is well suppressed.

Free-Falling Circular Cylinder in Quiescent Flow
The free-falling circular cylinder problem will be simulated in this subsection. In a 2 × 6 rectangular box, a circular cylinder located at (1,4) with a radius r = 0.125 is free falling from rest, as shown in Figure 18a. The densities of solid ( ) and fluid ( ) are 1.25 and 1, while the non-dimensional fluid viscosity and gravity for the investigated problem are 0 .1 and (0, −980). The cylinder positions and vorticity contours at different times are plotted in Figure 18. In order to validate the present framework, comparisons of the time history of cylinder positions and velocity are also plotted at Figure 19. It can be clearly shown that the preset results obtained by the proposed framework and methodology with ∆ = ∆ = 1/128 and ∆ = 10 −4 agree well with the reference benchmark results [29,30]. For the sake of completeness, CD, CL and are plotted in Figure 20a, while the effect of ve-

Free-Falling Circular Cylinder in Quiescent Flow
The free-falling circular cylinder problem will be simulated in this subsection. In a 2 × 6 rectangular box, a circular cylinder located at (1,4) with a radius r = 0.125 is free falling from rest, as shown in Figure 18a. The densities of solid (ρ s ) and fluid ρ f are 1.25 and 1, while the non-dimensional fluid viscosity and gravity for the investigated problem are 0 .1 and (0, −980). The cylinder positions and vorticity contours at different times are plotted in Figure 18. In order to validate the present framework, comparisons of the time history of cylinder positions and velocity are also plotted at Figure 19. It can be clearly shown that the preset results obtained by the proposed framework and methodology with ∆x = ∆y = 1/128 and ∆t = 10 −4 agree well with the reference benchmark results [29,30]. For the sake of completeness, C D , C L and ω are plotted in Figure 20a, while the effect of velocity modification method is plotted in Figure 20b. From Figure 20, it shows that the present framework can also well suppress for the free-falling cylinder problem. The pressure contours at t = 0.2 are also plotted to further show the benefit of employing velocity modification method. From Figure 21a, the smooth pressure is obtained, while Figure 21b clearly shows oscillated pressure near the downward cylinder interface.
1, while the non-dimensional fluid viscosity and gravity for the investigated problem are 0.1 and (0, −980). The cylinder positions and vorticity contours at different times are plotted in Figure 18. In order to validate the present framework, comparisons of the time history of cylinder positions and velocity are also plotted at Figure 19. It can be clearly shown that the preset results obtained by the proposed framework and methodology with ∆ = ∆ = 1/128 and ∆ = 10 −4 agree well with the reference benchmark results [29,30]. For the sake of completeness, CD, CL and are plotted in Figure 20a, while the effect of velocity modification method is plotted in Figure 20b. From Figure 20, it shows that the present framework can also well suppress for the free-falling cylinder problem. The pressure contours at t = 0.2 are also plotted to further show the benefit of employing velocity modification method. From Figure 21a, the smooth pressure is obtained, while Figure 21b clearly shows oscillated pressure near the downward cylinder interface.
To show the advantage of the present framework, different density ratios ( / ), 1.05 and 1.01, were also investigated in this study. Figure 22 clearly shows that the present framework can correctly and consistently predict solid velocity and position. The applicability and efficiency of the present framework for solving the free-falling cylinder problem is thus confirmed.      To show the advantage of the present framework, different density ratios (ρ s /ρ f ), 1.05 and 1.01, were also investigated in this study. Figure 22 clearly shows that the present framework can correctly and consistently predict solid velocity and position. The applicability and efficiency of the present framework for solving the free-falling cylinder problem is thus confirmed.

Drafting-Kissing-Tumbling (DKT) Problem of Two Free-Falling Circular Cylinders in Quiescent Flow
Finally, the drafting-kissing-tumbling (DKT) problem was investigated in this subsection. In a computational domain of ( , ) = (0~6, −1~1), there are two cylinders with a radius r = 0.125 located at (1, 0.001) and (1.5, −0.001). The density ratio and viscosity ratio was chosen as 1.5 and 0.01, while the gravity is (980,0). To model the DKT scenario in this study, the repulsion force model [29] was employed in this study, with the choice of a stiffness parameter as 5 × 10 −3 and effective force range as 3Δ . Due to the large repulsion force, the sub-time stepping method is employed to prevent the divergence issue when solving the equations of motions. The idea is to introduce a sub-time step Δ (= Δ / ), which is times smaller than the original time step Δ , and solves the sub-time step times. In practice, the fluid solver costs the most of the computational time, so the additional cost due to the sub-time stepping can be ignored. In this study, the sub-time step is set as 20. The predicted pressure and vorticity contours with Δ = Δ = 1/256 and Δ = 5 × 10 −5 at t = 0.15, 0.2, 0.25 and 0.3 are plotted in Figures 23 and 24. The comparisons of solid velocity and position with the referenced solution [2] are then made. From Figures 25 and  26, good agreements can be seen at ≤ 0.17, (drafting and first kissing stage). For > 0.17 (the later kissing and tumbling stage), it still shows reasonable agreements. As indicated by [2], the discrepancy should be due to strong instability and the choice of collision parameter and the initial horizontal offset of two cylinders. It is noted that there is no drafting for the present framework if the horizontal offset is set to zero.

Drafting-Kissing-Tumbling (DKT) Problem of Two Free-Falling Circular Cylinders in Quiescent Flow
Finally, the drafting-kissing-tumbling (DKT) problem was investigated in this subsection. In a computational domain of (x, y) = ( 0 ∼ 6, −1 ∼ 1), there are two cylinders with a radius r = 0.125 located at (1, 0.001) and (1.5, −0.001). The density ratio and viscosity ratio was chosen as 1.5 and 0.01, while the gravity is (980, 0). To model the DKT scenario in this study, the repulsion force model [29] was employed in this study, with the choice of a stiffness parameter as 5 × 10 −3 and effective force range as 3∆x. Due to the large repulsion force, the sub-time stepping method is employed to prevent the divergence issue when solving the equations of motions. The idea is to introduce a sub-time step ∆t s (= ∆t/s), which is s times smaller than the original time step ∆t, and solves the sub-time step s times. In practice, the fluid solver costs the most of the computational time, so the additional cost due to the sub-time stepping can be ignored. In this study, the sub-time step s is set as 20. The predicted pressure and vorticity contours with ∆x = ∆y = 1/256 and ∆t = 5 × 10 −5 at t = 0.15, 0.2, 0.25 and 0.3 are plotted in Figures 23 and 24. The comparisons of solid velocity and position with the referenced solution [2] are then made. From Figures 25 and 26, good agreements can be seen at t ≤ 0.17, (drafting and first kissing stage). For t > 0.17 (the later kissing and tumbling stage), it still shows reasonable agreements. As indicated by [2], the discrepancy should be due to strong instability and the choice of collision parameter and the initial horizontal offset of two cylinders. It is noted that there is no drafting for the present framework if the horizontal offset is set to zero.
additional cost due to the sub-time stepping can be ignored. In this study, the sub-time step is set as 20. The predicted pressure and vorticity contours with Δ = Δ = 1/256 and Δ = 5 × 10 −5 at t = 0.15, 0.2, 0.25 and 0.3 are plotted in Figures 23 and 24. The comparisons of solid velocity and position with the referenced solution [2] are then made. From Figures 25 and  26, good agreements can be seen at ≤ 0.17, (drafting and first kissing stage). For > 0.17 (the later kissing and tumbling stage), it still shows reasonable agreements. As indicated by [2], the discrepancy should be due to strong instability and the choice of collision parameter and the initial horizontal offset of two cylinders. It is noted that there is no drafting for the present framework if the horizontal offset is set to zero.

Conclusions
In this study, we proposed a framework to model the flow-particle interaction problem. The solid object is modeled by the DIIB method, while the fluid velocity and pressure are solved by the IDFC 2 method. To prevent the instability problem of a density ratio close to 1, the hydrodynamic force and torque which appeared in the equations of motion are directly evaluated with modified interpolation kernel function. The methodology of evaluating the drag/lift forces and the treatment of a circular cylinder approaching the wall have also been addressed. To validate the applicability and accuracy, problems of one and two particles were investigated. To make the framework more stable, the idea of a sub-time stepping method was also employed for DKT problems of two particles. From the investigating validation/benchmarking problems, the simulated results reveal the applicability, accuracy and efficiency of the present framework.