Consider coupled models for particulate flows, where the disperse phase is made of particles with distinct sizes. This leads to a system coupling the incompressible Navier–Stokes equations to the Vlasov–Fokker–Planck equations. For the model with random initial data near global equilibrium, we establish uniform regularity in suitable Sobolev spaces using energy estimates and demonstrate exponential decay of energy over 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 Stokes number, with exponential decay of an error over time, and propose a stochastic asymptotic-preserving (s-AP) scheme to simulate the behavior of multi-phase flow system, efficiently in both kinetic and hydrodynamic regimes. Additionally, we provide uniform error estimates of the bi-fidelity method for this coupled model with random initial inputs. Numerical examples illustrate the accuracy and efficiency of the method.