## 1 Introduction  
The past several decades have witnessed the rapid progress of the ability to simulate aerodynamic flows via computational fluid dynamics (CFD), which is believed to bring fundamental changes to the process of aerospace design [1]. However, the application of CFD in industry is still limited by the incapacity to simulate complex turbulent flows accurately. With the evolution of computing equipment, large eddy simulation (LES), though still largely confined in academia or research and development departments of industry, has proven to be increasingly useful in a wide range of applications in recent years such as aircraft engines flow, turbine combustor flame, and flows involving turbulent transition [2, 3, 4]. The LES approach stands between the Reynolds Averaged Navier-Stokes (RANS) and the Direct Numerical Simulation (DNS) methods in terms of accuracy and computational cost and has shown promising feasibility for vortex dominated flows found in various aerospace applications. In an LES computation, the large scale flow structures are directly simulated while the scales below a certain resolution are approximated by an explicit subgrid scale (SGS) model or represented implicitly by the inherent numerical dissipation. The explicit SGS models are commonly categorized into three types [5,6]: the functional models, which are designed to model the action of the subgrid terms and are expected to approximately dissipate energy in the smallest scales; the structural models, which are aimed at making the best approximation of the modeled terms and are focused on representing energy transfer between various length scales; and the mixed models, which are developed by combining a functional model and a structural model together as the latter alone is generally not dissipative enough. For the LES research a critical problem is to prevent subgrid scale models from being overwhelmed by discretization errors, if the numerical methods are not handled carefully. Previous investigation on finite volume method (FVM) concluded that the small scales of the flow might suffer from high numerical damping [7]. In terms of wave resolving, it was claimed that waves with wave lengths less than 10 times the filter width cannot be adequately resolved with a 2 nd order FVM because of the large truncation error of the mesh resolution, under the assumption that the filter width corresponds to the cell size [8]. To reduce the numerical error, plenty of high order methods such as discontinuous Galerkin (DG) [9], spectral difference (SD) [10] and flux reconstruction (FR) [11] or correction procedure via reconstruction (CPR) [12] methods have been developed in the last decades and have shown great potential in LES thanks to their lower numerical dissipation as well as their capacity of handling complex geometries. These methods have multiple degrees of freedom (DOFs) defined in their computing element according to the order of accuracy and are able to reduce truncation error if the order is set high enough, it is also possible for these methods to resolve shorter waves in the LES solution since their DOFs per wave (DPW) increase [8]. Plenty of successful and encouraging works including wall-modeled turbulence and transition have appeared in the literature [6,13, 14, 15, 16, 17, 18, 19, 20]. Among the aforementioned high order methods, the FR method provides a general high-order framework, which recovers the collocation-based nodal DG and SD methods in case of linear fluxes and also sets a path for new schemes. The FR approach is further developed into a family of energy-stable flux reconstruction (ESFR) schemes [21], whose stability has been proved for multiple types of elements [22, 23, 24]. It has been found that the FR approach is more computationally efficient than conventional lower order methods [25].