The presence of particles with a small but finite size, suspended in viscous fluids with low volumetric concentrations, is observed in many applications. The present study focuses on the tridimensional and incompressible lid-driven flow of Newtonian fluids through the application of the immersed boundary method and the Euler–Lagrange approach. These methods are used to numerically predict three-dimensional particle motion by considering nearly neutrally buoyant conditions as well as all relevant elementary processes (drag and lift forces, particle rotation, particle–wall interactions, and coupling between phases). Considering the current stage of the numerical platform, two coupling approaches between phases are considered: one-way and two-way coupling. A single particle is inserted in the cavity after steady-state conditions are achieved. Its three-dimensional motion is obtained from numerical simulations and compared with research data, considering the same conditions, evidently showing that the particle trajectory follows the experimental data until the first collision with a solid surface. After this first contact, there is a deviation between the results, with the two-way coupling results better representing the experimental data than the one-way coupling results. The dimensionless forces’ peaks acting on the particles are associated with the relative velocity of the particle near the wall–particle collision position. In terms of magnitude, in general, the drag force has shown greater influence on the particle’s motion, followed by the rotation-induced and shear-induced lift forces. Finally, a special application is presented, in which 4225 particles are released into the domain and their dynamic is evaluated throughout dimensionless time, showing similar behavior for both couplings between phases, with variations in local concentrations observed in certain regions. The mean square displacement used to quantify the dispersion evolution of the particles showed that the particulate flow reaches an approximately homogeneous distribution from the moment of dimensionless time tU/S = 130.