The numerical simulation of turbulent flows (for an introduction, see e.g. [Pope, 2000], [Tennekes and Lumley, 1972], [Mathieu and Scott, 2000], [Davidson, 2004]) is one of the great challenges in Computational Fluid Dynamics (CFD). It is commonly accepted that the physics of the flow of a continuous fluid is well represented by the Navier-Stokes equations. The direct numerical simulation (DNS) (for a review, see [Moin and Mahesh, 1998]) discretizes directly the three-dimensional Navier-Stokes equations. The basic requirement for such a simulation to succeed is the use of numerical schemes of high-order accuracy and meshes fine enough to capture the smallest scales of motion, to the order of the Kolmogorov scales. However, when the ratio of inertial forces to viscous ones, quantified by the Reynolds number, increases the smallest scales become smaller, and the amount of information (handled and processed) necessary for a Navier-Stokes based prediction becomes enormous. For homogeneous isotropic turbulence, for instance, the number of grid points needed is typically proportional to the 9/4 power of the Reynolds number. Then DNS applies well to rather academic problems characterized by simple geometry and low Reynolds numbers, and is a wonderful tool for the understanding of turbulence and validation of models. Nonetheless, the prediction of most engineering flows cannot be done by DNS with today’s computers, and probably not before computing power have been increased during a couple of decades.