In this work, high-fidelity large-eddy simulations are used to investigate hypersonic boundary layer transition. The initial motivation for this project stemmed from the joint research of CFDLAB Technion and Purdue University on hypersonic boundary layer transition over cone-cylinder-flare configurations. A blunt, axisymmetric, cone-cylinder-flare with 5-degree half-angle cone and 12-degree flare was proposed to perform experiments and numerical simulations on in hypersonic transitional flow. The compressible Navier-Stokes equations for a single, ideal gas were solved in generalized curvilinear coordinates. Careful consideration to preserve the Geometric Conservation Law and freestream preservation was taken. The convective terms of the governing equations were spatially discretized using the centralized Gradient-Based Reconstruction method to enable excellent turbulence resolution and sharp discontinuity capturing. The viscous terms of the governing equations were spatially discretized using an alpha-damping method to avoid well-known discretization errors. The third-order total-variation-diminishing Runge-Kutta method was used for time integration. To avoid infeasible grid point requirements at the Reynolds numbers considered in this study, the compressible equilibrium ODE wall-stress model was used to model the inner layer dynamics of the boundary layer. GPU acceleration was realized using the OpenACC programming standard. MPI was employed for parallelization amongst multiple GPUs. A suite of test cases was used to validate the implementation, compare against existing methods, and perform parametric studies. The solver was applied to Mach 6 transitional flow over the target geometry, yielding results that closely match experimental data in key aspects, including wall heat flux distributions, separation and reattachment locations, transition onset, and proper orthogonal decomposition (POD) modes.