CFD-DEM simulations have become a popular tool for scientists because they provide a means for accurately modeling particle migration in suspensions. These simulations use the Lagrangian framework in the DEM (Discrete Element Method) and the Eulerian framework in the CFD (Computational Fluid Dynamics) side. General CFD-DEM methods use a drag law approximation to calculate the force exerted on particles through the fluid. However, this method fails to accurately represent complex particle-fluid interactions in close distance between particles due to the lubrication effects. It is because the mesh resolution is not fine enough to capture the correct hydrodynamic interaction between particles at a distance less than 0.1 of particle diameter. As an alternative, a fully resolved simulation could capture the drag force on the particles without using any approximation due to their fine mesh resolution. [1-2] The present work used the same strategy by employing an Immersed Boundary Method (IBM) as a CFD solver in coupled CFD-DEM model. All recent progress in IBM was added to the code to increase its accuracy, especially in correctly capturing the lubrication force between particles. [3-5] Then, different test cases were considered to validate the code for short-range and long-range hydrodynamic interactions between particles. [6-7] Since IBM needs fine mesh resolution, its simulations take much longer than general CFD-DEM simulations. To overcome this matter, we modified the IBM code to implement it in parallel mode. By using parallel simulations, the behavior of multiple particles and their interactions were studied. Results confirmed that our approach could efficiently investigate the interactions between particles within the fluid phase to provide insight into the suspension characteristics.