## 1 Introduction  
The dynamics of the particle-wall and particle-particle collisions for Reynolds numbers below 1000 have been studied by various experiments and numerical simulations. [PERSON] et al. [1] and [PERSON] et al. [2] explored the vortex structures generated by particle-wall collisions by the experiment and numerical simulation, and [PERSON] et al. [3] analyzed the dynamics associated with particle-particle collisions. [PERSON] [4] and [PERSON] et al. [5] conducted numerical simulations of the flows around multiple particles to investigate the clustering of particles and turbulence effects. [PERSON] et al. [6], using the numerical simulation based on the Eulerian approach for the fluid phase and the Lagrangian approach for the solid phase, studied the effect of the diameter of dust particles injected in a spherical vessel on the resulting turbulent flow field and dust distribution. However, the flow structures and the particle behavior of a system with multiple particle-wall collisions are still unclear.  
The shot peening process, which involves impacting a metal surface with a large number of shots (particles), is widely used in practical mechanical engineering processes. When the particles collide with the metal surface, a compressible residual stress is generated over its surface [7]. The important parameter of this high-quality process is coverage level, which is defined as the percentage ratio of the total indentation area to the total treated surface area [8]. [PERSON] et al. [9] investigated the coverage level by noting the number of impressions. [PERSON] et al. [10] observed peer forming in experiments and numerical simulations. However, these studies focused on the plate, without adequately considering the interactions of the rebounding particles. [PERSON] et al. [11] investigated the behavior of particles in experiments and multi-body dynamics simulations, however, they largely ignored the influence of the flow structure. [PERSON] et al. [7] investigated the location of particle-surface impacts in the multiphase flow simulation. They employed the Euler-Lagrange approach and Reynolds-averaged Navier-Stokes equations, which cannot sufficiently model the particle-flow interactions and the detailed flow structures.  
Among the Cartesian mesh methods, the immersed boundary method (IBM) developed by [PERSON] et al. [12], has been widely used for its simplicity and applicability to multiple moving objects. We have developed an IBM-based two- and three-dimensional Euler-Euler compressible flow solver that analyzes the interaction between the particles and turbulence in a gas-particle flow [13; 14]. In the present study, we developed IBM-based three-dimensional incompressible flow solver [15] for exploring the structures and behaviors of particles undergoing multiple particle-wall collisions. This present flow solver directly performs a numerical simulation using the Euler-Euler approach. This simulation can capture the particle-flow interactions and the detailed flow structures.  
The aim of this paper is to investigate the accuracy of the present IBM flow solver and to examine the flows around multiple particle-wall collisions to elucidate the interaction between the flow and particles using the Euler-Euler approach. To investigate the accuracy of the IBM flow solver, we computed the flow structure and the drag coefficient around a fixed sphere using different schemes, and compared the results of the solver with those of existing methods. We analyzed the flow when the moving particles impinge on a wall. Finally, for examining the interaction between the flow and a lot of moving particles, we compared the two-way coupled simulation with the one-way simulation that is generally applied in clarifying the particle kinematics in industry and in investigating the flow around particle-wall collisions.