A 3D parallel high-order spectral difference (SD) solver with curved local mesh refinement is developed in this research to simulate flow through stenoses of varied degrees (50%, 60%, 65%, 70% and 75%) of radius constriction at inlet Reynolds number of 500. This solver employs high-order curved mesh in the vicinity of arterial wall and the local mesh refinement technique reduces the overall computational cost by distributing more elements in critical regions. In simulation of flow through stenosis of 50% radius constriction, velocity profiles predicted from the SD solver agree well with previous DNS results and experimental data. Mesh independency study shows that numerical results from a conforming and a non-conforming mesh agree well with each other. When the constriction degree is larger than 50%, visualizations through iso-surfaces of Q-criterion show that vortex rings are ejected from the stenosis throat, advecting downstream before they hit the vessel walls and they finally break down and merge into a large bulk region of small-scale turbulence. The observation is consistent with the vorticity contour which is characterized by development of the Kelvin-Helmholtz instability when shear layers are formed, rolled up and advected downstream between the central jet and the recirculation region. When the constriction degree turns to 75%, the flow transitions rapidly downstream of stenosis throat and dramatic pressure drop is witnessed. This provides a fluid-dynamic explanation for clinical definition of critical stenosis (i.e. over 75% luminal radius narrowing). Furthermore, pressure drop across a stenosis is found to be proportional to square of ratio of non-stenosed area to minimum area at the stenosis throat with a linear correlation coefficient equal to 0.9998. Finally, this solver is proven to have excellent scalability on massively parallel computers when multi-level refinement of meshes is performed to capture small-scale structures in the turbulence region.