Direct numerical simulations with a second-order immersed boundary-lattice Boltzmann method are used to investigate the effect of particle fluctuation on flows in fixed random arrays of spheres at low and moderate particle Reynolds numbers. Random velocities obeying the isotropic Maxwellian distribution are imposed to all particles in the computational domain to mimic the granular temperature due to particle collisions. The simulation results show that the effect of particle fluctuation on the drag force is significant especially when the solid volume fraction and the particle Reynolds number are small. The drag increment due to the particle fluctuation increases with the increase of the granular temperature-based Reynolds number, however, it decreases with the increase of the solid volume fraction and the particle Reynolds number. On the basis of the simulation results, a new drag force relation based on the stationary drag in random arrays of spheres for arbitrary solid volume fractions, granular temperature-based Reynolds numbers, and particle Reynolds numbers is formulated. The fluctuations of the drag force on individual particles with respect to the mean drag are also analyzed. It is found that the drag on individual particles can differ up to 40% from the mean at small granular temperature-based Reynolds numbers and the difference increases rapidly with the increase of the granular temperature-based Reynolds number, indicating the particle fluctuation can affect the individual drag force significantly. The present results also show the drag force on individual particles well follows Gaussian distribution. The mean relative deviation is found to increase with the increase of the granular temperature-based Reynolds number and the solid volume fraction but is not affected much by the particle Reynolds number. (C) 2017 Elsevier B.V. All rights reserved.