Computational fluid dynamics (CFD) is increasingly used to simulate complex industrial systems. Most CFD analysis relies on the Reynolds-averaged Navier-Stokes (RANS) approach and traditional two-equation turbulence models. Higher-fidelity approaches to the simulation of turbulence such as wall-resolved large eddy simulation (LES) and direct numerical simulation (DNS) remain limited to smaller applications or to large supercomputing platforms. Nonetheless, continued advances in supercomputing are enabling the simulation of physical systems of increasing size and complexity. These simulations can be used to gain unprecedented insight into the physics of turbulence in complex flows and will become more widespread as petascale architectures become more accessible. As the scale and size of LES and DNS simulations increase, however, the limitations of current algorithms become apparent. For larger systems, more temporal and spatial scale must be resolved, thus increasing the time-scale separation. While the smaller time scales dictate the size and the computational cost associated with each time step, the larger time scales dictate the length of the transient. An increased time-scale separation leads to smaller time steps and longer transients, eventually leading to simulations that are impractical or infeasible. In practice the presence of multiple and strongly separated time scales limits the effectiveness of CFD algorithms for LES and DNS applied to large industrial systems. Moreover, the situation is likely to become worse as even larger systems are simulated, thus increasing the size and length of transients. At the same time transients currently simulated on petascale architectures are unlikely to become any faster on exascale architectures. This paper presents an ensemble-averaging technique for transient simulations, aimed at collecting averaged turbulent statistics faster. The focus is on ergodic flows and simulations. Ensemble averaging involves creating multiple models and combining them to produce a desired output. This technique is commonplace in machine learning and artificial neural networks, and it is at the basis of RANS/URANS turbulence modeling. In the proposed approach, multiple instances of the same ergodic flows are simulated in parallel for a short time and summed to create an ensemble. Provided each instance is sufficiently statistically decorrelated, this allows considerable reduction in the time to solution. This paper focuses on the theory and implementation of the methodology in Nek5000, a massively parallel open-source spectral-element code. Also presented is the application of the method to the DNS and LES simulation of channel flow and pipe flow.