Consider coupled models for particulate flows, where the disperse phase is made of particles with distinct sizes. We are thus led to a system coupling the incompressible Navier–Stokes equations to the Vlasov–Fokker–Planck equations. For the model with random initial data near the global equilibrium, a uniform regularity is established in some suitable Sobolev spaces by using energy estimates and we prove that the energy decays exponentially in time by hypocoercivity arguments. We also prove that the generalized polynomial chaos stochastic Galerkin (gPC-sG) method has spectral accuracy, uniformly in time and the Knudsen number, and the error decays exponentially in time. A numerical scheme is designed to simulate the behavior of the multi-phase system. This scheme is asymptotic-preserving, which allows the efficient computation in both kinetic and hydrodynamic regimes. Numerical examples illustrate the accuracy and asymptotic behavior of the scheme, with several interesting applications.