The objective of this study is to develop direct numerical simulations (DNS) to investigate the aerodynamic performance, transition to turbulence, and to capture the laminar separation bubble occurring on a wind turbine blade. Simulations are conducted with spectral/hp element method to investigate the details of flow separation bubble over wind turbine blades with NACA-4412 airfoil at wide range of design parameters. This airfoil is chosen because recent studies have shown that it is challenging to capture the details of the flow instabilities and pressure fluctuations in the separated shear layer of wind turbines by experimental methods. Furthermore, owing to more accurate development of DNS, the separated bubbles at high Reynolds numbers are captured. The results show that the vortex structures shed from the trailing edge of the airfoil by raising the angle of attack (α). Consequently, the fully turbulent flow develops downstream of the trailing edge (Karman vortex). Moreover, the pressure fluctuation significantly increased by raising α. However, some rolling up of the flow structures, similar to Kelvin–Helmholtz rolls, on the pressure surface near the trailing edge, are observed at α>12 °. The separation point was delayed from X sep/C = 0.19 to 0.58 by decreasing α from 16 to 0 at Re = 5 × 10 4.