Computational analyses of fluid flow through packed pebble bed domains using the Reynolds-averaged Navier–Stokes (RANS) framework have had limited success in the past. Because of a lack of high-fidelity experimental or computational data, optimization of Reynolds-averaged closure models for these geometries has not been extensively developed. In the present study, direct numerical simulation (DNS) was employed to develop a high-fidelity database that can be used for optimizing Reynolds-averaged closure models for pebble bed flows. A face-centered cubic (FCC) domain with periodic boundaries was used. Flow was simulated at a Reynolds number of 9308 and cross-verified by using available quasi-DNS data. During the simulations, low-frequency instability modes were observed that affected the stationary solution. These instabilities were investigated by using the method of proper orthogonal decomposition, and a correlation was found between the time-dependent asymmetry of the averaged velocity profile data and the behavior of the highest-energy eigenmodes. Finally, the effects of the domain size and the method of averaging were investigated to determine how these parameters influenced the stationary solution. A violation of the ergodicity assumption was observed.