This paper presents a new cosmological hydrodynamical simulation method based on a moving mesh defined by the Voronoi tessellation of a set of discrete points. The method uses a finite-volume approach to solve the hyperbolic conservation laws of ideal hydrodynamics with a second-order unsplit Godunov scheme and an exact Riemann solver. The mesh-generating points can move arbitrarily, and when stationary, the scheme is equivalent to an ordinary Eulerian method with second-order accuracy. When moving with the local flow, it becomes a Lagrangian formulation that is fully Galilean-invariant, unlike traditional Eulerian codes. This property is crucial for cosmological simulations involving supersonic bulk flows. The method automatically adjusts spatial resolution and inherits the principal advantage of SPH for simulating cosmological structure growth. It retains the high accuracy of Eulerian methods in shock treatment and improves contact discontinuity handling. The method is implemented in the new code AREPO, which is parallelized for distributed memory computers and allows adaptive refinement or derefinement of the unstructured mesh. It also includes an individual time-step approach for finite volume hydrodynamics and a high-accuracy treatment of self-gravity for gas, enabling seamless coupling with high-resolution collisionless dark matter simulations. The paper discusses the performance of AREPO through a suite of test problems, demonstrating its effectiveness as a competitive alternative to current SPH and Eulerian techniques. The method is based on Delaunay and Voronoi tessellations, which are efficient for mesh generation and allow for continuous mesh deformation without mesh-tangling effects. The paper also discusses the algorithms for generating these tessellations, including incremental insertion, projection of convex hulls, recursive subdivision, and flipping. The method is robust and efficient, with a data structure that allows rapid access to topological objects of the tessellation. The paper also addresses the treatment of degenerate cases, where points lie on edges or faces of the current tessellation, requiring special handling to ensure correct geometric predicates. The method uses exact long integer arithmetic to evaluate geometric tests when necessary, ensuring robustness and accuracy. The paper concludes that the moving-mesh approach represents a compromise between SPH and AMR, combining the automatic adaptivity and geometric flexibility of SPH with the high-accuracy treatment of shocks, shear waves, and fluid instabilities of AMR. The new code AREPO is shown to be a promising and competitive method for future applications in cosmology and other fields.This paper presents a new cosmological hydrodynamical simulation method based on a moving mesh defined by the Voronoi tessellation of a set of discrete points. The method uses a finite-volume approach to solve the hyperbolic conservation laws of ideal hydrodynamics with a second-order unsplit Godunov scheme and an exact Riemann solver. The mesh-generating points can move arbitrarily, and when stationary, the scheme is equivalent to an ordinary Eulerian method with second-order accuracy. When moving with the local flow, it becomes a Lagrangian formulation that is fully Galilean-invariant, unlike traditional Eulerian codes. This property is crucial for cosmological simulations involving supersonic bulk flows. The method automatically adjusts spatial resolution and inherits the principal advantage of SPH for simulating cosmological structure growth. It retains the high accuracy of Eulerian methods in shock treatment and improves contact discontinuity handling. The method is implemented in the new code AREPO, which is parallelized for distributed memory computers and allows adaptive refinement or derefinement of the unstructured mesh. It also includes an individual time-step approach for finite volume hydrodynamics and a high-accuracy treatment of self-gravity for gas, enabling seamless coupling with high-resolution collisionless dark matter simulations. The paper discusses the performance of AREPO through a suite of test problems, demonstrating its effectiveness as a competitive alternative to current SPH and Eulerian techniques. The method is based on Delaunay and Voronoi tessellations, which are efficient for mesh generation and allow for continuous mesh deformation without mesh-tangling effects. The paper also discusses the algorithms for generating these tessellations, including incremental insertion, projection of convex hulls, recursive subdivision, and flipping. The method is robust and efficient, with a data structure that allows rapid access to topological objects of the tessellation. The paper also addresses the treatment of degenerate cases, where points lie on edges or faces of the current tessellation, requiring special handling to ensure correct geometric predicates. The method uses exact long integer arithmetic to evaluate geometric tests when necessary, ensuring robustness and accuracy. The paper concludes that the moving-mesh approach represents a compromise between SPH and AMR, combining the automatic adaptivity and geometric flexibility of SPH with the high-accuracy treatment of shocks, shear waves, and fluid instabilities of AMR. The new code AREPO is shown to be a promising and competitive method for future applications in cosmology and other fields.