High-order numerical methods have proven to be an essential tool to improve the accuracy of simulations involving turbulent flows through the solution of conservation laws. Such flows appear in a wide variety of industrial applications and its correct prediction is crucial to reduce the power consumption and improve the efficiency of these processes. The present study implements and analyzes different types of high-order spatial discretization schemes for unstructured grids to assess and quantify their accuracy in simulations of turbulent flows. In particular, high-order Finite Volume methods (FVM) based on least squares and fully constrained deconvolution operators are considered and their accuracy is evaluated in a variety of linear and non-linear test cases and throughanalytical analysis. Special emphasis is placed on the comparison of formally second-order and high-order FVM, showing that the former can over-perform the latter in terms of accuracy and computational performance in under-resolved configurations. High-order Spectral Element methods (SEM), including Spectral Difference (SD) and Flux Reconstruction (FR), are compared in different linear and non-linear configurations. Furthermore, a SD GPU-based solver (based on the open-source PyFR solver) is developed and its performance with respect to other state of the art CPU-based solvers will be discussed, showing that the developed GPU-based solver outperforms other state of the art CPU-based solvers in terms of performance-per-euro and performance-per-watt. The accuracy and behavior of SEM under aliasing are assessed in linear test cases using analytical tools. The use of grids with high-order cells, which allow to better describe the surfaces of interests of a given simulation, in combination with SEM is also analyzed. The latter analysis demonstrates that special care must be taken to ensure appropriate numerical accuracy when utilizing meshes with such elements. This document also presents the development and the analysis of the Spectral Difference Raviart-Thomas (SDRT) method for two and three-dimensional tensor product and simplex elements. This method is equivalent to the SD formulation for tensor product elements and it can be considered as a natural extension of the SD formulation for simplex elements. Additionally, a new family of FR methods, which is equivalent to the SDRT method under certain circumstances, is described. All these developments were implemented in the open-source PyFR solver and are compatible with CPU and GPU architectures. In the context of high-order simulations of turbulent flows found in rotor-stator interaction test cases, a sliding mesh method (which involves non-conformal grids and mesh motion) specifically tailored for massivelyparallel simulations is implemented within a CPU-based solver. The developed method is compatible with second-order and high-order FVM and SEM. Grid movement, needed to simulate rotor-stator test cases due to the relative movement of each domain zone, is treated using the Arbitrary-Lagrangian-Eulerian (ALE) formulation. The analysis of such formulation depicts its important influence on the numerical accuracy and stability of numerical simulations with mesh motion. Moreover, specific non-conformal discretization methodscompatible with second-order and high-order FVM and SEM are developed and their accuracy is assessed on different non-linear test cases. The parallel scalability of the method is assessed with up to 11000 cores, proving appropriate computational efficiency. The accuracy of the implementation is assessed through a set of linear and non-linear test cases. Preliminary results of the turbulent flow around a DGEN 380 fan stage in an under-resolved configuration are shown and compared to available experimental data.