Abstract : The simulation of high-Reynolds number massively separated flows is a challenging problem of prime interest in industry. Even in this context we have to investigate numerical tools which are able to work in a Large Eddy Simulation mode as well as RANS or the hybrid RANS-LES mode. It is then necessary to use an adapted numerical technology: in the present work, the spatial discretization is based on a mixed finite element/finite volume formulation on unstructured grids. The numerical dissipation of the upwind scheme is made of sixth-order space derivatives and corresponds to a fifth-order term in trunca-tion error, in order to limit as far as possible the interactions between numerical and subgrid scale (SGS) dissipation, . The turbulence model is based on a hybridization strategy which blends a variational mul-tiscale large-eddy simulation (VMS-LES) equipped with dynamic SGS models and a two-equation RANS model. Particular attention is paid to the VMS-LES approach used in this work. The separation between the large and the small resolved scales is obtained through a variational projection operator based on spatial average on agglomerated cells . The dynamic procedure (Germano) allows the adaptation of the constant of the SGS model to the spatial and temporal variation of the flow characteristics, while the VMS formulation restricts the SGS model effects to the smallest resolved scales . The hybridization strategy uses a blending parameter, such that a VMS-LES simulation is applied in region where the grid resolution is sufficiently fine to resolve a significant part of the turbulence fluctuations, while a RANS model is used in the regions of coarse grid resolution . For instance, near the body surface, it is used under two forms, either as a statistical model for the turbulent boundary layer or as a (thinner) wall law for the LES part of the flow. The capability of the proposed hybrid model to accurately predict the aerodynamic forces acting on a circular cylinder in the supercritical regime and tandem cylinders is investigated. It is shown that rather coarse grids can be used to obtain a reasonable prediction of important bulk coefficients.