  Hydrodynamic cosmological simulations at present usually employ either the
Lagrangian SPH technique, or Eulerian hydrodynamics on a Cartesian mesh with
adaptive mesh refinement. Both of these methods have disadvantages that
negatively impact their accuracy in certain situations. We here propose a novel
scheme which largely eliminates these weaknesses. It is based on a moving
unstructured mesh defined by the Voronoi tessellation of a set of discrete
points. The mesh is used to solve the hyperbolic conservation laws of ideal
hydrodynamics with a finite volume approach, based on a second-order unsplit
Godunov scheme with an exact Riemann solver. The mesh-generating points can in
principle be moved arbitrarily. If they are chosen to be stationary, the scheme
is equivalent to an ordinary Eulerian method with second order accuracy. If
they instead move with the velocity of the local flow, one obtains a Lagrangian
formulation of hydrodynamics that does not suffer from the mesh distortion
limitations inherent in other mesh-based Lagrangian schemes. In this mode, our
new method is fully Galilean-invariant, unlike ordinary Eulerian codes, a
property that is of significant importance for cosmological simulations. In
addition, the new scheme can adjust its spatial resolution automatically and
continuously, and hence inherits the principal advantage of SPH for simulations
of cosmological structure growth. The high accuracy of Eulerian methods in the
treatment of shocks is retained, while the treatment of contact discontinuities
improves. We discuss how this approach is implemented in our new parallel code
AREPO, both in 2D and 3D. We use a suite of test problems to examine the
performance of the new code and argue that it provides an attractive and
competitive alternative to current SPH and Eulerian techniques. (abridged)
