The efficient utilization of oil shale resource is beneficial to environment and the resource diversity. In this work, a bubbling fluidized bed (BFB) reactor was first used to pyrolyze Huadian oil shale experimentally, mainly focusing on the influence of fluidization flow rates. The yield of shale oil, char and non-condensable gas varies little at different fluidization flow rates, but the components in shale oil are different due to the secondary cracking, such as the fractions and species of aliphatic hydrocarbons, aromatic hydrocarbons and non-hydrocarbons. Then, the dense discrete phase model (DDPM) was used to track the trajectories of oil shale particles with different diameters in a bubbling fluidized bed. Population balance model (PBM) coupled Eulerian model was used to well describe the bubbling status of bed materials to characterize the hydrodynamic behaviors. Furthermore, a segmented kinetic model was compiled by UDF module to illuminate the pyrolysis of oil shale particles. The multi-phase reactive flow modeling clearly indicated the distribution of oil shale particles and gaseous products in BFB reactor, as well as pyrolysis time and the secondary cracking at different fluidization flow rates. In summary, this work provides a reference basis for complete pyrolysis of oil shale particles in a bubbling fluidized bed and generating shale oil with high yield and good quality.