This paper is devoted to numerical modeling of a supersonic flow around a blunt body by a viscous gas with an admixture of relatively large high-inertia particles that, after reflection from the surface, may go beyond the shock layer and change the flow structure dramatically. To calculate the gas-dynamic interaction of moving particles with the shock layer, it is important to take into account the large difference in scales of the flow around the particles and around the body. To make the computations effective, we use a meshless method to solve non-stationary Navier–Stokes equations. The algorithm is based on the approximation of partial derivatives by the least squares method on a set of nodes distributed in the calculation area. Each moving particle is surrounded by a cloud of calculation nodes belonging to its domain and moving with it in space. The algorithm has been tested on the problem of the motion of a single particle and a pair of particles in a supersonic flow around a sphere.