Understanding the velocity field is very important for modern cosmology: it gives insights to structure formation in general, and also its properties are crucial ingredients in modelling redshift-space distortions and in interpreting measurements of the kinetic Sunyaev-Zeldovich effect. Unfortunately, characterising the velocity field in cosmological N-body simulations is inherently complicated by two facts: i) The velocity field becomes manifestly multi-valued after shell-crossing and has discontinuities at caustics. This is due to the collisionless nature of dark matter. ii) N-body simulations sample the velocity field only at a set of discrete locations, with poor resolution in low-density regions. In this paper, we discuss how the associated problems can be circumvented by using a phase-space interpolation technique. This method provides extremely accurate estimates of the cosmic velocity fields and its derivatives, which can be properly defined without the need of the arbitrary "coarse-graining" procedure commonly used. We explore in detail the configuration-space properties of the cosmic velocity field on very large scales and in the highly nonlinear regime. In particular, we characterise the divergence and curl of the velocity field, present their one-point statistics, analyse the Fourier-space properties and provide fitting formulae for the velocity divergence bias relative to the non-linear matter power spectrum. We furthermore contrast some of the interesting differences in the velocity fields of warm and cold dark matter models. We anticipate that the high-precision measurements carried out here will help to understand in detail the dynamics of dark matter and the structures it forms.

17 pages, 17 figures, matches version accepted by MNRAS; minor changes: improved discussion of real space properties and relation to smoothed fields, explicit derivation of multi-stream derivatives