A simulation method is proposed to explore the effect of particle size polydispersity in magnetorheology including Brownian motion. The method aims to extend the classical particle-level simulation methodology developed by Klingenberg et al. [J. Chem. Phys. 91, 7888-7895 (1989)] for the case of polydisperse magnetorheological (MR) fluids. The simulation study concerns the aggregation kinetics at rest as well as the rheological behavior under start-up of steady shear and dynamic oscillatory shear tests at increasing strain amplitudes. Results demonstrate that the effect of polydispersity is only relevant at the transition regime between magnetostatic to hydrodynamic control of the suspension structure. The yielding behavior is correlated to the structural characteristics (radial distribution functions, pair correlation functions, and angular connectivities) of the MR fluids before the onset of flow. A more abrupt transition is observed for polydisperse MR fluids because interparticle links are weaker in this case if compared to monodisperse suspensions in spite of the fact that polydisperse MR fluids exhibit a larger connectivity. (C) 2015 The Society of Rheology.